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ABSTRACT 

Motivated by recent suggestions that many Be stars form through binary mass transfer, we searched the APOGEE survey for Be 
stars with bloated, stripped companions. From a well-defined parent sample of 297 Be stars, we identified one mass-transfer binary, 
HD 15124. The object consists of a main-sequence Be star (Mge = 5.3 + 0.6 Mo) with a low-mass (Maonor = 0.92 + 0.22 Mo), 
subgiant companion on a 5.47-day orbit. The emission lines originate in an accretion disk caused by ongoing mass transfer, 
not from a decretion disk as in classical Be stars. Both stars have surface abundances bearing imprint of CNO processing in 
the donor's core: the surface helium fraction is Yye ~ 0.6, and the nitrogen-to-carbon ratio is 1000 times the solar value. The 
system's properties are well-matched by binary evolution models in which mass transfer begins while a 3 — 5 Mo donor leaves 
the main sequence, with the originally less massive component becoming the Be star. These models predict that the system will 
soon become a detached Be - stripped star binary like HR 6819 and LB-1, with the stripped donor eventually contracting to 
become a core helium-burning sdO/B star. Discovery of one object in this short-lived (^1 Myr) evolutionary phase implies the 
existence of many more that have already passed through it and are now Be + sdO/B binaries. We infer that (10 — 60) % of Be 
stars have stripped companions, most of which are ~ 100 x fainter than the Be stars in the optical. Together with the dearth of 
main-sequence companions to Be stars and recent discovery of numerous Be + sdO/B binaries in the UV, our results imply that 
binarity plays an important role in the formation of Be stars. 
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1 INTRODUCTION synthesis models that include binaries make significantly different 
predictions from single-star models for the ionizing flux, metal en- 
richment, and supernova rates of unresolved stellar populations (e.g. 
Stanway et al. 2016; Ma et al. 2016; Steidel et al. 2016), the under- 
lying binary evolution models have been subject to only limited tests 
with observations of individual binaries in the nearby Universe (e.g. 


Han et al. 2003; Sana et al. 2013; Sen et al. 2021). 


Binary mass transfer underlies many of the most important open 
questions in astrophysics. A large fraction of all stars — particularly 
massive stars — are in binaries that exchange mass at some point in 
their evolution (e.g. Sana et al. 2012; Kobulnicky et al. 2014; Moe & 
Di Stefano 2017), and mass transfer often dramatically changes the 
lives and deaths of stars (e.g., de Mink et al. 2013). Realistic mod- 


eling of interacting binaries is one of the major areas of research in 
modeling of unresolved stellar populations (e.g. Eldridge et al. 2017). 
Mass transfer is challenging to study observationally because it often 
Occurs on short timescales, with most observed binaries found before 
or after periods of mass transfer. Thus, although stellar population 
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A somewhat counterintuitive but common outcome of binary evo- 
lution with mass transfer is mass ratio inversion (e.g. Giuricin et al. 
1983), wherein the less-massive component of a binary is observed to 
be the more evolved star. Such “Algol-type” binaries form when the 
initially more massive component of a binary loses a large fraction 
of its mass to a companion (e.g. Paczyński 1971), generally through 
stable (i.e., long-lived and non-dynamical) mass transfer. Explaining 
the existence of Algol-type binaries as a consequence of mass trans- 
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fer was one of the great achievements of stellar evolution modeling 
in the last century (see Pustylnik 1998, for a review). 

Although stably mass-transferring binaries have been studied in- 
tensely for the better part of a century, several recent developments 
have renewed interest in their evolution. The first is the mounting evi- 
dence that spin-up through mass transfer is among the primary causes 
of the “Be” phenomenon; i.e., the prevalence of rapidly-rotating B 
stars with emission lines, which usually trace circumstellar disks. 
The idea that Be stars might be binary evolution products spun up by 
mass transfer, first popularized by Pols et al. (1991), has recently been 
bolstered by the discovery in the UV of more than a dozen, compact 
core-helium burning stars (*sdO/B" stars) as companions to Be stars 
(e.g. Wang et al. 2018, 2021). Some of these companions have also 
been confirmed interferometrically (Klement et al. 20212). Other in- 
direct lines of evidence have also recently suggested that many Be 
stars have hard-to-detect companions stripped by binary interactions: 
the frequency of main-sequence companions to Be stars is unusually 
low (Bodensteiner et al. 20202), and the SEDs of many Be star disks 
hint at tidal truncation by a companion (Klement et al. 2019). Several 
lines of evidence also suggest that sdO/B stars are essentially always 
formed via binary interactions (e.g. Han et al. 2003; Pelisoli et al. 
2020; Geier et al. 2022), so studying the formation channel of Be 
stars can also shed light onto the formation of sdO/B stars. 

Second, several unusual binaries have recently been identified that 
most likely contain a Be star and a warm (Teg ~ 15,000 K), bloated 
(R « 5 Ro) companion (Liu et al. 2019; Rivinius et al. 2020; Saracino 
et al. 2021). All of these systems were initially interpreted as con- 
taining a main-sequence B star orbiting a black hole, but subsequent 
studies (e.g. Irrgang et al. 2020; Shenar et al. 2020; Bodensteiner 
et al. 2020b; El-Badry & Quataert 2021; El-Badry & Burdge 2021) 
have shown that the B stars are not on the main sequence, but appear 
to be undermassive (~ 0.5— 1.5 Mo) stripped products of binary evo- 
lution that are currently contracting to become core helium burning 
sdO/B stars. Because the current evolutionary state of these objects is 
short-lived, their discovery implies that there may be many other sys- 
tems with similar evolutionary histories yet to be discovered. These 
Be + sdO/B binaries would be long-lived but challenging to identify: 
the Be star usually contributes the vast majority of the total light in 
the optical, and the predicted few-km s^! radial velocity (RV) shifts 
of the Be star are only marginally detectable due to Be stars’ large 
rotation velocities and disk-driven spectral variability. 

The chances of detecting compact Be + sdO/B systems are higher 
in the UV, where the hot sdO/B star contributes a larger fraction of 
the total light. Even so, only the tail of the population with the hottest 
sdO/B stars are unambiguously detectable (e.g. Schootemeijer et al. 
2018). Compared to typical sdO/B stars, a majority of the sdO/B 
companions to Be stars discovered in the UV (e.g. Gies et al. 1998; 
Koubsky et al. 2012; Peters et al. 2013; Wang et al. 2018, 2021) 
are unusually hot and luminous. This is at least partly a selection 
effect, since companions are more likely to be detectable when they 
contribute a larger fraction of the total light. 

While sdO/B companions are most easily detected in the UV, the 
optical and IR wavelengths are optimal for characterizing earlier 
evolutionary stages of systems that will later become Be + sdO 
binaries. For instance, shortly after leaving the main sequence, the 
initially more massive star can be relatively luminous and rich in 
spectral features in the optical and/or IR regions, and all but invisible 
in the UV. Surveys in different wavelength regimes are thus necessary 
to build a more complete picture of the evolutionary paths that can 
lead to rapidly rotating Be stars and their stripped companions. 

In this paper, we further test the hypothesis that binary mass trans- 
fer is an important formation channel for Be stars by conducting a 
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systematic search for Be stars in binaries being spun up by mass 
transfer. This allows us to then obtain a more robust estimate of the 
fraction of Be stars that have gone through a mass-transfer phase than 
could be obtained from HR 6819 and LB-1, as both systems were 
identified serendipitously from ill-defined effective parent samples. 
The search yields one binary, HD 15124, on which most of the paper 
is focused. 

The remainder of this paper is organized as follows. Section 2 
describes our search for Be + stripped star binaries in APOGEE. 
Section 3 presents data on HD 15124, including the APOGEE spec- 
tra, follow-up optical spectra, light curves, radial velocities, and the 
broadband SED. Section 4 summarizes our constraints on the physi- 
cal parameters of the binary. In Section 5, we construct binary evo- 
lution models to investigate the system’s formation history and pre- 
dicted future evolution. We discuss the implications of our results 
for the broader Be star population in Section 6. The paper’s main 
results are summarized in Section 7, and the Appendices provide 
more information on several aspects of our modeling. 


2 APOGEE SEARCH FOR BE STARS WITH 
NARROW-LINED COMPANIONS 


We began by searching DR16 of the APOGEE survey (Majewski et al. 
2017; Jonsson et al. 2020) for Be stars with narrow-lined companions. 
Although the primary goals of the APOGEE survey are centered on 
cool giants, it also observed several thousand OB stars all across the 
sky, mainly for use as telluric standards. 


2.1 Identifying B and Be stars 


As a first step, we selected hot and luminous stars observed 
by APOGEE in the Gaia color-magnitude diagram (CMD). Of 
the 437,445 stars included in DR16, 432,541 had a correspond- 
ing source (as identified in the allStar table) in Gaia DR2. 
339,374 of these had at least a marginally significant parallax 
(parallax_over_error > 3), and 338,733 also had a reported 
bp_rp color. From this sample, we selected candidate OB stars as 
those satisfying Mg < 0.5, where Mg = G + 5log (w/100), and 
bp_rp < 1 (gray points in Figure 1). This left us with 6,176 hot 
and luminous star candidates with APOGEE spectra. 5,295 of these 
were observed at APO, 1,018 were observed at LCO, and 359 were 
observed at both sites. 

Colored lines in Figure 1 show synthetic photometry for single-star 
MIST evolutionary models (Choi et al. 2016) with Ay = 0 (solid) 
and Ay = 1 (dot-dashed). These show that the sample is expected 
to contain sources with masses M z 3 Mo and moderate extinction. 
For example, a typical 5 Mo main-sequence star with Teg = 16000 K 
and log g = 4 would have MG ~ —1 and Gpp — Ggp ~ —0.2 with no 
extinction, and would thus fall within our selection region as long as 
Ay < 1.75. This CMD-based selection is crude, and more careful 
selections exist in the literature (e.g. Zari et al. 2021). Our goal was 
to cast a wide net and include a majority of hot and luminous stars 
in our sample, removing cooler stars later. 

We retrieved the combined spectra (apStar/asStar files) of all 
6,176 hot and luminous star candidates and inspected them indi- 
vidually to identify Be stars. We focused on the Brackett 114 
line at 16,802 A, which is usually the strongest emission line in Be 
stars within the APOGEE wavelength coverage (Chojnowski et al. 
2015). This search yielded 297 Be stars with unambiguous emission 
in the Brackett lines; these are shown in Figure | with gold stars. 
For comparison, we also show the sample of 261 Be star candidates 
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Figure 1. Sample of candidate Be stars in context. Upper left: Color-magnitude diagram. Gray points show APOGEE DR16 sources with Ggp — Grp < | and 
Mg < 0.5; i.e., hot and luminous stars. Colored lines show MIST evolutionary tracks without extinction (solid) and with Ay = 1 (dot-dashed). Upper right: 
We inspected the APOGEE spectra of the 6176 sources and identified 297 Be stars with significant emission in the Brackett series (gold stars). Open stars show 
another sample of Be stars identified from APOGEE spectra by Chojnowski et al. (2015), a majority of which overlap with our sample and a few of which fall 
outside our CMD selection region. Red star shows HD 15124. Lower left: sky distribution. Most Be stars, including HD 15124, are in the Galactic plane. Lower 
right: apparent magnitude and distance. HD 15124 is in the nearest and brightest quartile of the Be star sample. 


identified from APOGEE DR12 by Chojnowski et al. (2015) as open 
star symbols. Examples of several spectra are shown in Figure 2. 


146 of the Be stars in our sample are also in the Chojnowski et al. 
(2015) sample. A majority of our new identifications are stars ob- 
served after DR12. Of the 115 Be stars identified by Chojnowski 
et al. (2015) that are not included in our sample, 18 are not in- 
cluded in the DR16 allStar table, and 16 did not pass our Gaia 
parallax_over_error cut. 36 passed other cuts, but fell outside 
our CMD selection (left panel of Figure 1). Finally, 45 fell within 
our selection region in the CMD, but were not identified as Be stars 
in our visual inspection. We re-inspected the spectra of these 45 
stars and found that most of them show weak or ambiguous emission 
(e.g., bottom panel of Figure 2) and therefore were not flagged in our 
search. 


As is evident from the upper right panel of Figure 1, 14% of the Be 
stars from the Chojnowski et al. (2015) sample fall outside our CMD 
selection region. Most of these sources have high extinctions. This 
suggest that we could have obtained a 10-20% larger sample of Be 


stars by expanding our CMD selection region. We opted for a more 
restrictive selection because a larger CMD search region would lead 
to more contamination from giants and lower-mass main sequence 
stars, significantly increasing the number of spectra requiring visual 
inspection. 


In summary, the parent sample within which we search for bloated 
stripped companions contains 297 Be stars. The sample is bright 
(median apparent magnitude G = 9.4) and relatively nearby (me- 
dian distance 1.3 kpc). Most stars have precise parallaxes (median 
parallax_over_error = 33 with Gaia eDR3 data). After correct- 
ing for extinction and reddening using the 3D dust map from Green 
et al. (2019), the median color and absolute magnitude of the sample 
is (bp rp)g = —0.07 and Mg o = —1.6. This corresponds to spectral 
type B3; i.e., ~ 5 Mo for stars near the main sequence. This is a 
slightly lower mass than the spectral type where the frequency of 
Be stars peaks (B1.5-B2; e.g. Kogure & Hirata 1982; Rivinius et al. 
2013), reflecting the fact that lower-mass stars are longer-lived and 
favored by the IMF. 
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22 Search for Be stars with narrow lines and/or RV variability 


For each of the 297 Be stars, we retrieved and inspected the 
individual-visit spectra (apVisit/asVisit files) for evidence of 
binarity. The median number of visits per target in DR16 is 5. All 
but a few targets have at least 3 visits, and about a third have more 
than 7. The median time baseline spanned by the observations for a 
single target is 295 days, and the median single-visit SNR is 208 per 
pixel. 
The primary features we searched for were: 


(1) Epoch-to-epoch RV variability: For low-mass companions with 
periods < 1 year, as expected for binaries with ongoing or recently 
terminated mass transfer, we expect the companions to have typical 
observable epoch-to-epoch RV shifts of > 10km s^!. The RV shifts 
of the Be star, on the other hand, are generally not expected to be 
detectable, due to the Be stars' large masses and rapid rotation rates. 
We found that the RVs reported in the allVisit table are often unre- 
liable for Be stars, so we did not use these in our analysis. Instead, we 
searched for objects in which there appeared to be a genuine Doppler 
shift of photospheric absorption lines after barycentric correction. 

(ii) Narrow absorption lines: Be stars are expected to be rapidly 
rotating, and thus, to have rotationally smeared-out absorption lines 
except when viewed nearly pole-on. In contrast, the bloated, stripped 
companions we seek to identify are expected to have been tidally 
synchronized by mass transfer, and thus to have relatively low pro- 
jected rotation velocities. We thus regarded the presence of narrow 
absorption lines as a promising indicator of a binary companion. 

(ili) Line profile variability: We also searched for epoch-to-epoch 
changes in emission line profile shapes, which could also be related 
to binarity. Such changes are, however, ubiquitous in the emission 
lines of Be stars (e.g. Slettebak 1979; Dachs et al. 1981; Okazaki 
1991) and are thought to be due primarily to changes in the structure 
of the circumstellar disk, often unrelated to binarity. We therefore 
focused primarily on absorption lines, and did not consider emission 
line profile changes alone to be strong evidence for binarity. 


2.2.1 Search sensitivity 


The initial motivation for our search was to identify relatively hot, 
bloated companions similar to those discovered in LB-1 and HR 
6819. The stripped stars in these objects have already terminated 
mass transfer and begun contraction toward the extreme horizontal 
branch, with effective temperatures of 13,000 and 16,000 K. 

Unfortunately, stars with Teg z; 8000 K rapidly lose their narrow 
metal lines in the near infrared (Figure B1), making such companions 
very challenging to detect with APOGEE. One might still hope to 
detect the Brackett line cores of an RV-variable companion shifting 
from epoch to epoch, but this is not easily realized in practice because 
the emission line profiles of Be stars are themselves highly variable. 

As we show in Appendix B, our search is thus sensitive to com- 
panions with T.g < 8000K. Fortunately, the stripped companions 
in systems like LB-1 and HR 6819 are expected to have been cooler 
earlier in their evolution, when mass transfer was still ongoing. Emis- 
sion lines are still expected to exist during this stage, but they will 
originate in an accretion disk, not in a decretion disk as in classical 
Be stars. This is the evolutionary stage to which the search is most 
sensitive. 


2.2.2 Search results 


Our visual inspection resulted in three objects being flagged as po- 
tential Be star binaries. 
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The first was HD 38616, whose spectra contain narrow lines but 
showed no evidence of RV variability between two visits separated by 
3 days (3rd panel of Figure 2) . We obtained several high-resolution 
follow-up spectra with the FEROS spectrograph at La Silla (Kaufer 
et al. 1999) over a 7 day period in December 2020. These confirmed 
the object to be a Be star with narrow absorption lines (spectral type 
B3) but ruled out RV variability of the absorption lines with ampli- 
tude larger than 1 km s7! over the 7-day baseline. We conclude that 
this object is likely a normal "shell star"; i.e., a Be star viewed nearly 
edge-on, with significant self-absorption from the circumstellar disk 
(e.g. Struve 1931; Rivinius et al. 2006). As one of the few Be stars in 
the TESS continuous viewing zone, it remains an interesting object 
to study in detail. 

The 2nd object of interest we identified was HD 55606, whose 
spectra showed significant epoch-to-epoch line profiles variations in 
many emission lines and unusually strong metallic emission lines, but 
no obvious narrow absorption lines (4th panel of Figure 2). It turned 
out that this object was already studied in detail by Chojnowski et al. 
(2018) and Wang et al. (2021), who found it to be a Be + sdO binary; 
i.e., a Be + stripped star binary in which the donor star has already 
contracted to become a hot subdwarf. 

The final candidate was HD 15124, on which most of our anal- 
ysis is focused. We describe the data for the system in Section 3 
and constrain physical parameters in Section 4. Our constraints are 
summarized in Table 1. 


3 DATA FOR HD 15124 
3.1 Previous literature 


HD 15124 was classified as spectral type B5 III (Tag ~ 15000 K) by 
Davis et al. (1973) and as B3/4 V/IV (Teg = 16000 — 18000 K) by 
Jensen (1981). An R « 3,000 optical spectrum of the object was 
obtained by Huang et al. (2010), who inferred Tag = 18603 + 300K 
and logg = 4.21 + 0.04 by fitting the Balmer lines, and vsini = 
87 + 13 from the He I lines. The APOGEE ASPCAP pipeline found 
very different spectral parameters, T.g = 6200 + 65 K and logg = 
3.76 + 0.05. None of these analyses accounted for the presence of 
a second star or disk, so we suspect their reported uncertainties are 
significantly underestimated. Chojnowski et al. (2015) classified the 
object as a Be star based on its APOGEE spectra. To our knowledge, 
no previous studies recognized the object as a spectroscopic binary. 


3.2 APOGEE spectra 


HD 15124 was observed 11 times by the APOGEE sur- 
vey between October 2016 and November 2017 (APOGEE ID 
2M02281038+5716293). The spectra have a typical SNR of 500 
per pixel, wavelength coverage of 15150 to 16950 (with the standard 
two chip gaps), and spectral resolution R ~ 22,500. Figure 3 shows 
three 250 A-wide cutouts of the APOGEE spectra from three visits, 
which are chosen to span the range of RV shifts in the data. 

In contrast to the vast majority of Be stars, HD 15124 has nar- 
row absorption lines, suggesting the presence of a slowly rotating 
star (v sini ~ 20kms7!; Section 4.5). Almost all the lines in the 
APOGEE spectra undergo coherent RV shifts from epoch to epoch, 
with a maximum shift of ~ 145kms~!. The exceptions are (a) dif- 
fuse interstellar bands (DIBs) at 15615, 15651 and 15671 A, and (b) 
double-peaked emission in the hydrogen Brackett lines (Br 11, 12, 13, 
14, and 15 at 16811, 16411, 16114, 15885, and 15705 À, with emis- 
sion becoming weaker in the higher-order lines). As we will show, 


A Be star is born 5 


T T 
BD 4- 622345 : normal B star 


= 
O 


normalized flux 
e 
Co 
i 


HD 221188: normal Be star 


0.8 J 


normalized flux 


T 
HD 38616: shell star 


1.0 k 


0.8 J 


normalized flux 


1.4F mp55606: Be--sdO « |l ] 


PS 
: | 
m | 
g2 | | 
5 
T 1 4 
BE 1.0 — WS hU 
o i Li 
a V 
0.8 -3 | | | | 
v BD + 62108: Chojnowski + 15 Be star; not included in this work 
3 
= 1.0 ial "T" ! 
"2 
D 
= 
E: 
B 08 | 
o 
S 


16550 16600 16650 16700 16750 16800 16850 16900 
wavelength [A] 


Figure 2. APOGEE spectral cutouts of representative B stars examined during our search. The strongest feature is the Br11 line at 16811 A. Different colors 
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Figure 3. APOGEE spectra of HD 15124. We show spectra from three different visits, which are labeled by orbital phase (Section 3.7). All narrow absorption 
lines shift coherently, with a ~140kms~! shift between the visits plotted in blue and red. In addition to RV-variable absorption lines, double-peaked emission 
is evident in the hydrogen Brackett series (most obvious in Br11 at 16811 À). RV shifts of the emission lines are weaker than those of absorption lines. Several 
strong lines are labeled; these all trace the narrow-lined, RV-variable donor star. 


HD 15124 contains a cool, slowly-rotating star that is losing mass to of continuum dilution by the Be star, which contributes more than 
a more-massive, hotter, rapidly rotating star with an accretion disk. half the flux in the H-band, but no narrow lines. 

We refer to the slowly-rotating, obviously RV-variable component as 

the “donor”, and the accreting companion as the “Be star". 


The flux calibration of APOGEE visit spectra is not expected to 3.3 Low-resolution, flux calibrated Shane/Kast spectrum 


be precise (Nidever et al. 2015), so most of our analysis is based We observed HD 15124 using the Kast spectrograph (Miller & Stone 


on pseudo-continuum normalized spectra. We define the pseudo- 1994) on the 3m Shane telescope at Lick Observatory on November 3, 
continuum using a spline fit to regions of the spectra without obvious 2020. We used the 600/7500 grating on the red side and the 600/4310 
absorption lines. This does not necessarily correspond to the “true” grism on the blue side, with the D55 dichroic and a 1 arcsec slit. 
continuum (which is in any case ill-defined for a system containing This resulted in wavelength coverage of 3300-8400 À with spectral 
two RV-variable stars and a disk), but simply allows us to bring resolution R ~ 1, 000. We observed the A3V standard star HD 37725 
observed spectra and models to a common scale. for flux calibration. A 120-second exposure yielded SNR « 300 per 
The narrow metal lines in the APOGEE spectra are relatively pixel at 4000 À. We reduced the spectrum using pypeit (Prochaska 
shallow, with the strongest lines having normalized depths ranging et al. 2020). 
from 0.05 to 0.1. As we will show, this is in large part a consequence The spectrum is shown in Figure 4. Most narrow metal lines from 
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Figure 4. Flux-calibrated optical spectrum (R « 1000) of HD 15124 and spectral models. Top panel shows our best-fit model spectra, with the contributions of 
the Be star (blue) and donor (orange) shown separately. The Be star dominates the optical spectrum, with the donor contributing ~ 20%. Bottom panels show 
the effects of extinction (left; we constrain Ay based on the Green et al. 2019 3D dust map), the temperature of the Be star (middle; it is constrained by the 
strength of the Balmer jump, and the Be star radius is adjusted to match the normalization at 4000 A for each temperature), and helium abundance (right; the 


strong observed He I lines require significant He enhancement). 


the donor are not detectable in the spectrum due to the low spectral 
resolution; most of the obvious features come from the Be star. A 
notable exception is the Ca II H&K lines, which are strong in the 
donor but weak in the hotter Be star. Although its spectral resolution is 
low, the Kast spectrum provides strong constraints on the temperature 
and radius of the Be star. In particular, the Balmer jump at 3646 À is 
sensitive primarily to the Be star's effective temperature, with weak 
dependence on its surface gravity and composition. With only a low- 
resolution spectrum, it would be difficult to determine how much 
of the optical flux comes from the donor. Fortunately, the donor's 
contributions to the spectrum are constrained reliably from the high- 
resolution spectrum. 


3.4 High-resolution LBT/PEPSI spectrum 


We observed HD 15124 for 300 seconds on September 28, 2020 us- 
ing the Potsdam Echelle Polarimetric and Spectroscopic Instrument 
(PEPSI; Strassmeier et al. 2015) spectrograph on the Large Binocular 


Telescope in binocular mode. The spectrum was reduced as described 
in Strassmeier et al. (2018); it covers the wavelength ranges of 4755- 
5428 A and 6236-7433 A with spectral resolution R ~ 50,000 and 
per-pixel SNR ~ 500. The spectrum was obtained at phase ¢ = 0.59 
(see Section 3.7). Figure 5 shows portions of the spectrum and our 
best-fit model (Section 4). Zoomed-in cutouts of the spectrum are 
also shown in Figures 10 and 13. 

Like the APOGEE spectra, the PEPSI spectrum contains double- 
peaked emission lines. Emission is obvious in Ha and more subtle, 
but still unambiguously present, in HG. There is also double-peaked 
emission in some, but not all, of the He I lines, with a typical sep- 
aration of ~ 300 km s^! between the peaks. There are no detectable 
Fe lines in emission. 

The PEPSI spectrum makes it obvious that at least two luminous 
stars contribute in the optical, because it contains features charac- 
teristic of both hot and cool stars. In particular, the presence of 
many narrow metal lines (all due to the donor) suggests a temper- 
ature near 7000 K, while the presence of moderately strong helium 
line in absorption (due to the Be star) suggests a temperature of at 
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Figure 5. PEPSI spectrum of HD15124 (black) and best-fit binary model (red; see Section 4 for details). Greyed-out regions of the spectra are masked due to 
contamination from tellurics and interstellar absorption. Bottom three panels show regions of the spectra with prominent emission lines from the disk. Most of 
the visually obvious absorption lines are from the cool donor star (orange line in the upper panel of Figure 4), which contributes only ~ 20% of the light. The 


contributions of the Be star are most evident in the Balmer and He I lines. 


least 15,000 K. Fitting the PEPSI and Kast spectra simultaneously 
(Section 4) allows us to constrain the atmospheric parameters and 
abundances of both components. 


3.5 NRES spectra 


We observed HD 15124 with the NRES spectrograph (Siverd et al. 
2018) on the 1 mtelescope at Wise observatory through the Las Cum- 
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bres Observatory Global Telescope network (Brown et al. 2013) 32 
times between February 1 and March 8 of 2021. The wavelength 
coverage is 3800-8200 A, with spectral resolution R ~ 30, 000. Ex- 
posure times ranged from 600 to 1800 seconds, yielding a typical 
per-pixel SNR of 5 at 5,000 Å. The spectra were reduced with the 
Banzai-NRES pipeline!. The low SNR of these spectra made them 


1 https://github.com/LCOGT/banzai-nres 


Table 1. Physical parameters and 1o uncertainties for both components of 


HD 15124. 


Parameters of the unresolved source 


Right ascension o [deg] 37.043254 
Declination 6 [deg] 57.274803 
Apparent Gaia eDR3 magnitude G [mag] 7.98 
Apparent Gaia eDR3 color Gpp — Grp [mag] 0.42 
Corrected Gaia eDR3 parallax w [mas] 1.64 + 0.03 
Gaia ruwe 1.11 

Color excess E(g — r) [mag] 0.255 x 0.02 
Dereddened absolute magnitude Me, [mag] —1.62 + 0.06 
Parameters of the subgiant donor star 

Effective temperature Teg [K] 6900 + 250 
Surface gravity log(g/(cms7?)) 2.8 +0.15 
Projected rotation velocity v sini [kms-!] 20 +3 
Rotation velocity vrot [km s^] 51 x9 
Synchronization parameter Vrot/ (2x R/ Py) 0.94 x 0.17 
Microturbulent velocity Vmic [km s^! ] 2.5 £1 
Cont. flux ratio at 5000 A füonor/ ftot 0.20 + 0.02 
Cont. flux ratio at 16000 A fonol fat 0.35 + 0.03 
Radius R [Ro] 5.82 + 0.45 
Bolometric luminosity log(L/Lo) 1.83 + 0.09 
Mass M [Mo] 0.92 + 0.22 
Parameters of the Be star 

Effective temperature Teg [K] 16000 + 1000 


Surface gravity log(g/(cms"?)) 4.1 € 0.2 
Projected rotation velocity v sini [kms-!] 95 x 10 

Cont. flux ratio at 5000 À foel fiot 0.76 + 0.02 
Cont. flux ratio at 16000 Å Joel fiot 0.45 + 0.05 
Radius R [Ro] 4.1 +0.2 

Mass M [Mo] 5.3 + 0.6 
Bolometric luminosity log( L/ Lo) 3.00 x 0.12 
Rotation velocity Vrot 241 + 30 
Fraction of critical rotation Vrot/ Verit 0.60 + 0.08 
Parameters of the binary 

Orbital period Porb [day] 5.4692 + 0.0001 
Conjunction time To [MHJD] 57678.740 + 0.001 
Donor RV semi-amplitude Kaonor [km s^] 74.3 x 0.9 
Center-of-mass velocity yg [km s^! ] -11.1 x 0.7 
Donor eccentricity Edonor < 0.02 

Donor RV scatter Sdonor [km s7!] « 1.5 

Orbital inclination i [deg] 23.2 x 1.5 
Separation a [Ro] 24.0 + 0.8 

Mass function f (Mpe) [Mo] 0.23 + 0.01 


Surface abundances — assumed equal for both components 


Helium mass fraction Yue 0.63 + 0.08 
Carbon (relative to solar) C -2.3 + 0.2 
Nitrogen (relative to solar) N 0.8 + 0.3 
Oxygen (relative to solar) O —0.8 + 0.2 
Tron (relative to solar) Fe 0.05 + 0.1 
Be star disk parameters 

Cont. flux ratio at 5000 A Jäisk/ fiot 0.04 + 0.02 
Cont. flux ratio at 16000 A Jäisk/ fiot 0.2 x 0.05 
Emission line peak separation AVpeak [km s] 300 + 50 
Outer disk radius Router [Ro | 7.0 x 3.5 
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useful primarily for measuring RVs, which we used to constrain the 
binary's orbit (Section 3.7 and Figure 8). 


3.6 Light curves: Hipparcos, WISE, TESS, and KELT 


To search for photometric variability, we retrieved light curves of HD 
15124 from several time-domain photometric surveys. Unfortunately, 
the star is so bright that its photometry is saturated in many time- 
domain surveys (e.g. ASAS-SN, ZTF, CSS). However, the available 
photometry still allows us to understand the object's photometric 
variability over a range of timescales and wavelengths. 


3.6.1 Hipparcos 


The Hipparcos mission (ESA 1997) produced 132 usable photomet- 
ric measurements of HD 15124 (HIP 11487) over a 3 year baseline 
with median uncertainty of 0.013 mag. We retrieved the light curve 
from the ESA archive (Disk 2, Volume 17) and show it in the upper 
left panel of Figure 6. The sparse sampling of the data make it difficult 
to probe short-timescale variability, but longer-timescale variability 
is obvious. Most notably, the source brightened by about 0.05 mag 
in the first half of 1991, and then faded again by ~0.02 mag in 1992. 


3.6.2 WISE 


HD 15124 was also observed by the WISE mission (Wright et al. 
2010). We retrieved multi-epoch photometry from the AIIWISE Mul- 
tiepoch Photometry database (Cutri et al. 2021), which includes us- 
able light curves in both the W1 and W2 bands. The W2-band light 
curve, which has typical uncertainty 0.022 mag, is shown in the upper 
right panel of Figure 6. The behavior of the W1 light curve is similar, 
but its uncertainties are larger. Like the Hipparcos light curve, the 
WISE data shows clear evidence of long-term flux variability, with 
the mean magnitude brightening by 0.09 mag over a half-year period. 

Long-timescale brightness variability in both the optical and in- 
frared is ubiquitous in Be stars (e.g. Rivinius et al. 2013). In most 
cases, variability is attributed to growth and shrinking of a circum- 
stellar disk. Variability is typically more pronounced in the infrared, 
because the disk contributes a larger fraction of the total light there. 
The amplitude of the optical and NIR variability observed for HD 
15124 (x 0.1 mag) is on the low side for classical Be stars, where 
variability amplitudes of ~ 0.5 mag are common (e.g. Labadie-Bartz 
et al. 2017). The accretion disk in HD 15124 — being truncated by 
a close companion — must be smaller than the disks in classical Be 
stars, which often extend to tens of stellar radii. 


3.6.3 TESS 


HD 15124 was observed by the TESS satellite (Ricker et al. 2015) 
for 25 days during sector 18 in November 2019. We used eleanor 
(Feinstein et al. 2019) to perform background subtraction and PSF 
photometry on the TESS full-frame images, extracting a light curve 
with 30 minute cadence. We experimented with a variety of apertures, 
finding the shape of the extracted light curve relatively robust as long 
as the central pixels nearest to HD 15124 are included. The star 
has no comparably bright neighbors within several arcminutes, so 
contamination from background sources is expected to be modest. 
The extracted light curve, which spans 25 days and includes a 
3-day gap for data downlinking, is shown in the middle left panel 
of Figure 6. We calculate the phase of the TESS observations using 
the orbital ephemeris from the RVs (Section 3.7), which is defined 
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Figure 6. Light curves for HD 15124. Top panels show data from Hipparcos (over 3 years) and WISE (two sets of observations separated by half a year). On 
long timescales, the source's brightness varies at the 0.05-0.1 mag level, likely due to changes in the structure of the Be star's accretion disk. Center left panel 
shows the TESS light curve, which spans 25 days. The data is phased based on the RV ephemeris (Section 3.7). We overplot a model light curve calculated with 
PHOEBE for the stellar parameters calculated in Section 4 (not a fit). The variability in the PHOEBE light curve is primarily due to reflection and ellipsoidal 
variation, as illustrated schematically in the center right panel. Bottom panels show the mesh corresponding to the PHOEBE light curve at three different phases. 
Reflection (i.e., irradiation of one side of the donor by the Be star) dominates the short-timescale variability. 


such that the Be star would eclipse the donor (if the inclination were 
close to edge-on) at phase 0. The TESS light curve reveals complex 
variability on short timescales. Most pronounced is quasi-sinusoidal 
variability on the orbital period, with semi-amplitude of (1-2)%. Ad- 
ditionally, there is short-timescale (hours to days) variability that 
is manifest as jagged, somewhat irregular fluctuations in the light 
curve, as well as longer-timescale variations. We do not find com- 
parable variability in the TESS light curves of most other nearby 
bright stars, suggesting that both the short- and long-timescale vari- 
ability observed in HD 15124 is astrophysical. Irregular photometric 
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variability is common in accretion disks due to a combination of 
turbulence and changes in disk structure (e.g. Bruch 1992). 


3.6.4 KELT 


HD 15124 was observed by the Kilodegree Extremely Little Tele- 
scope (KELT; Pepper et al. 2007) between 2012 and 2015, yielding 
1355 150-second exposures. The data were obtained using a 42 mm- 
aperture telescope at Winer Observatory with a 26 x 26 degree field 
of view and a broad photometric filter that transmits most light red- 
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Figure 7. KELT optical light curve of HD 15124. Left panel shows the raw light curve, which contains ~1300 measurements over 2.5 years. A few points 
with much fainter magnitudes, likely due to systematics, fall outside the axis limits. Middle panel shows a zoom-in on a 5-month period (30 orbits), revealing 
variability on a range of timescales. Right panel shows the full light curve, phased to the ephemeris from the RVs. The phased light curve has similar shape to 
the TESS light curve (Figure 6), revealing a combination of ellipsoidal variability and the reflection effect. 


ward of 5000 À. Along with raw light curves, the KELT pipeline 
produces a detrended light curve, which filters out longer-timescale 
variability and systematics common with other nearby stars using the 
Trend Filtering Algorithm (TFA; Kovács et al. 2005). HD 15124 falls 
in KELT field N17, the photometry of which is as yet unpublished. 

The KELT light curve is shown in Figure 7. To approximately con- 
vert the KELT instrumental magnitudes to Johnson-Cousins R-band, 
we subtract 4.1 from the reported values. Like the TESS data, the 
raw light curve (left and middle panels) reveals both short- and long- 
timescale variability. To reduce scatter due long-term variability, we 
crudely detrended the raw light curve using a running median filter 
with width of 33 days (6 orbital periods).? We also clipped outliers 
with calibrated KELT magnitudes outside the range (7.90, 8.05). The 
resulting phase-folded light curve (right panel) has a similar shape 
to the TESS light curve, suggesting that the light curve shape is 
relatively stable over long timescales. 


3.6.5 Origin of the photometric variability 


Also plotted in Figure 6 is a model light curve computed with 
PHOEBE (Prša & Zwitter 2005; Prša et al. 2016; Horvat et al. 2018), 
a program which builds on the Wilson & Devinney (1971) code for 
modeling binary light curves. The model we show is not a fit to the 
observed light curve, but is a prediction based on the binary param- 
eters estimated from the spectra and RVs (Section 4). Although it 
does not reproduce all the observed variability, the model does ex- 
plain the most prominent feature in the light curve: variability on the 
orbital period, with minima that are sharper than maxima, almost 
resembling eclipses. 

This light curve shape results from the combined effects of reflec- 
tion/irradiation and tidal distortion of the donor. The irradiated, Be 
star-facing side of the donor is hotter and brighter than the opposite 
side. This alone leads to approximately sinusoidal variability on the 
orbital period. In addition, the donor star is tidally distorted into a 


? We also experimented with using the TFA detrended light curve. This 
yielded similar results, with somewhat less scatter, but also with a ~20% 
lower variability amplitude, suggesting that the detrending removed some 
real astrophysical variability. 


"teardrop" shape, such that the geometric cross-section visible to an 
observer changes over the course of the orbit. This also produces 
nearly sinusoidal variability, but with a period of half the orbital pe- 
riod. Adding these two effects together leads to the light curve shape 
predicted by PHOEBE, with sharp minima and flat maxima (middle 
right panel of Figure 6)? The origin of this variability is illustrated 
further in the bottom panels of Figure 6, which show the scene used 
by PHOEBE to calculate the model light curve at different phases. 
Due to reflection and ellipsoidal distortion, the surface temperature 
of the donor star is predicted to vary significantly (by about 1000 K) 
over the star's surface, with the hottest point being the side facing the 
Be star. This means that the flux-weighted mean temperature (i.e., 
what we measure spectroscopically) is expected to vary somewhat 
with orbital phase. This is indeed what we find from the APOGEE 
spectra: the spectroscopic effective temperature is ~ 500 K hotter at 
phase 0 than at phase 0.5. 

As we discuss in Section 4.4, the inclination of HD 15124 is 
close to face-on (i ~ 23%). The amplitude of variability due to 
reflection and ellipsoidal variability scale respectively as sini and 
sin? i (e.g. Faigler & Mazeh 2011), so the observed variability is 
significantly weaker than would be expected for a similar binary 
viewed close to edge-on. The additional light from the Be star, which 
is not significantly irradiated or tidally distorted, further dilutes the 
variability amplitude. 


3.7 Radial velocities 


We measured radial velocities of the donor in HD 15124 by cross- 
correlating a model spectrum with both the APOGEE spectra and 
the NRES spectra. We only consider the NRES spectra from which 
we obtained RV uncertainties of less than 20 km sl, as the rest were 
obtained in poor conditions and are consistent with noise. This left us 
with 11 RVs from APOGEE and 29 from NRES, spanning a 4 year 


5 The light curve morphology — both the periodic and irregular behavior — 
are similar to that of the mass-transfer binary HD 63021 recently studied by 
Whelan et al. (2021). We suspect that the origin of the photometric variability 
in that system, which Whelan et al. did not explain, is also a combination of 
reflection and ellipsoidal variability. 
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baseline. We were unable to measure RVs for the Be star, because 
its contributions to the APOGEE spectra are subtle (primarily excess 
emission in the Brackett series), and the NRES spectra do not have 
sufficient SNR to extract reliable RVs from the helium absorption 
lines. 

The resulting RVs are listed in Table A1 and shown in the upper 
panels of Figure 8. We fit them with a Keplerian orbit using a com- 
bined simulated annealing Markov chain Monte Carlo method, as 
described in El-Badry et al. (2018b). We also fit an RV "scatter" term 
C'scatter. Which we add in quadrature to the measured RV uncertainties 
to represent either underestimated uncertainties or intrinsic scatter 
in the RVs. Given the large number of RVs, the posterior is reason- 
ably well-behaved and the period is unambiguously determined to 
be Porp ~ 5.47 days. We first left the eccentricity free and found 
e & 0.02 + 0.01; i.e., marginal evidence of eccentricity. Due to the 
short orbital period and tidal deformation of the donor, we expect 
the orbit to be tidally circularized and suspect that the nonzero ec- 
centricity favored by the initial fit is simply a manifestation of RV 
scatter. We therefore fixed the eccentricity to 0, leading to an in- 
ferred Oscatter % 1 km s7 l| Nonzero scatter (with a larger magnitude) 
was also found in HR 6819 (Bodensteiner et al. 2020b; El-Badry & 
Quataert 2021), where it is likely due to non-radial pulsations in the 
donor. The scatter in HD15124 might have similar origin, but might 
also stem from long-term RV variability due a tertiary (Appendix C), 
or differences between the NRES and APOGEE RV zeropoints. 

The orbital solution of the donor yields a joint constraint on the 
component masses and inclination, as described in Section 4.4. 


3.8 Broadband SED 


To estimate the contributions of the circumstellar disk as a func- 
tion of wavelength, we compared the source's measured fluxes from 
broadband photometry to synthetic photometry predicted based on 
the parameters of the donor and Be star alone. Johnson-Cousins 
U BV photometry for HD 15124 was measured from the ground by 
Deutschman et al. (1976), while near-IR and IR photometry was 
obtained by the 2MASS (Skrutskie et al. 2006) and WISE (Wright 
et al. 2010) missions. No UV observations of HD 15124 are pub- 
lished from GALEX (Morrissey et al. 2007). Fortuitously, a UV flux 
was measured by one of the earliest UV missions, the Celescope 
experiment (Davis 1968). We converted the measured Celescope 
magnitude, U3 = 11.53 mag (Davis et al. 1973) to a flux density 
using the empirical calibration from Houziaux (1974). This yielded 
a well-sampled SED from 0.1 to 20 um, which is plotted in Figure 9. 
We inflate all the reported photometric uncertainties by adding a 
10% flux uncertainty in quadrature to the reported values to con- 
servatively account for HD 15124's photometric variability and the 
unknown phasing of the photometry. 

We predict bandpass-averaged mean magnitudes for both com- 
ponents using empirically-calibrated theoretical models from the 
BaSeL library (v2.2; Lejeune et al. 1997, 1998). We assume a Cardelli 
et al. (1989) extinction law with total-to-selective extinction ratio 
Ry = 3.1, and we use an extinction prior from the Green et al. 
(2019) 3D dust map. We use pystellibs* to interpolate between 
model spectra, and pyphot? to calculate synthetic photometry. In 
the bottom panel of Figure 9, we plot the ratio of the observed fluxes 
and those predicted from the model of the donor and Be star, neglect- 
ing the disk. The observed UV and UBV fluxes are consistent with 


^ http;//mfouesneau.github.io/docs/pystellibs/ 
5 https://mfouesneau.github.io/docs/pyphot/ 
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coming only from the Be star and donor. This is almost by construc- 
tion, since we required the model spectra to match the flux calibrated 
blue-optical spectrum (Figure 4). However, there is clear evidence of 
excess emission at NIR and IR wavelengths, as is very common for 
Be stars. 


3.9 Astrometry 


After correcting for the parallax zeropoint (Lindegren et al. 2021) and 
underestimated parallax uncertainties for bright stars (El-Badry et al. 
2021c), the Gaia eDR3 parallax of HD 15124 is w = 1.64 x 0.03 
(eDR3 source id 458048851856088064; Gaia Collaboration et al. 
2021). This corresponds to a distance d = (610 + 11) pc. The re- 
ported ruwe value is 1.11, indicating a reasonably good single-star 
astrometric solution. For our best-fit parameters (Section 4), the pro- 
jected semimajor axis of the binary is p = © x (a/1 au) ~ 0.17 mas, 
about 1196 of the parallax. This may account for the somewhat larger 
than average ruwe (e.g. Belokurov et al. 2020) but is small enough 
that we expect the parallax to be reliable. 

Although the Gaia astrometric solution appears to be well- 
behaved, comparison of the Hipparcos-Gaia position change with 
the measured Gaia proper motion reveals evidence of long-term ac- 
celeration of the photocenter (Kervella et al. 2021). We discuss this 
in detail in Appendix C, where we conclude that HD 15124 may bea 
triple with a faint tertiary. The expected mass of the putative tertiary 
(^ 1 Mo) is low enough that it would be expected to contribute less 
than one percent of the total luminosity, so it has little effect on our 
modeling of the system. 


4 PARAMETERS OF THE BINARY 


To constrain the atmospheric parameters, radii, and abundances of 
both stars, as well as the distance and extinction, we simultane- 
ously fit the flux-calibrated low-resolution spectrum (Figure 4) and 
the normalized high-resolution spectra (Figures 3 and 5). The flux- 
calibrated spectrum primarily constrains the temperature and radius 
of the Be star, while the high-resolution spectrum constrains the 
donor atmospheric parameters, radii, projected rotation velocities, 
and abundances of both stars. 


4.1 Spectral fitting and uncertainties 


We simultaneously fit the Kast, PEPSI, and APOGEE spectra. We 
generated an irregular grid of 1D LTE atmospheres and synthetic 
spectra using ATLAS 12 (Kurucz 1970, 1979, 1992) to compute the 
atmosphere structure and SYNTHE (Kurucz 1993) for the radiative 
transfer calculations. We re-computed the atmosphere structure for 
each combination of atmospheric parameters and abundances. We 
use the linelist maintained by R. Kurucz.Ó To interpolate between 
models and predict spectra with atmospheric parameters and abun- 
dances not represented in the grid, we used a 2nd-order polynomial 
spectral model (e.g. Rix et al. 2016; Ting et al. 2019). Models were 
generated at resolution R = 300, 000, broadened using the rotation 
kernel from Gray (1992), and convolved with a Gaussian line spread 
function with spectral resolution (FWHM) R = 50,000 for PEPSI 
spectra, R = 22,500 for APOGEE spectra, and R = 1,000 for Kast 
spectra. 


6 http://kurucz.harvard.edu/linelists.html 
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Figure 8. Upper panels: radial velocities of the donor. Black and cyan point show RVs measured from APOGEE spectra and follow-up observations with NRES. 
Right panel shows phased RVs and the corresponding orbital solution. Bottom panels show the dynamically-implied mass of the Be star, given the observed 


mass function, f (Mpe) = 0.23 Mo, for a range of donor masses and inclinations. Gray shaded region shows plausible Be star masses, given the observed 
temperature and radius (Figure 11). Together, these constrain the orbital inclination to i ~ 23.2 + 1.5 degrees. Hatched region shows the donor mass constraint 


if the donor fills its Roche lobe (Figure 12). 


We pseudo-continuum normalize the APOGEE and PEPSI spec- 
tra using a 3rd order spline fit to regions without strong absorption 
lines, and we normalize the model spectra in the same way. When 
predicting model normalized spectra for the unresolved binary, we 
add unnormalized model fluxes for the two stars prior to normaliza- 
tion, as described in El-Badry et al. (2018a). For the model APOGEE 
spectra, we include 20% continuum dilution from the Be star’s disk, 
as suggested by the SED (Figure 9). We do not normalize the Kast 
spectrum, which sets the flux scale, and thus the radii of both stars. 


We began with a spectral model with 10 labels: the effective tem- 
peratures, surface gravities, radii, and projected rotation velocities 
of both stars, the microturbulent velocity vturb of the donor, and a 
global metallicty [Fe/H], varying all metals in lockstep with [Fe/H] 
and assuming the solar abundance pattern. It quickly became clear 
that this was not a sufficient set of abundance labels, because (a) 
the model helium lines were weaker than observed, and (b) carbon 
lines in the observed spectra are extremely weak, while other metal 
lines are relatively strong. This suggested that the abundances are 
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Figure 9. Broadband spectral energy distribution. Points with error bars show measurements from the Celescope UV mission (“U3”; Davis et al. 1973), 
ground-based UBV photometry (Deutschman et al. 1976), 2MASS (JH Ks; Skrutskie et al. 2006), and WISE (W 1W2W 3W 4; Wright et al. 2010). Blue 
and red lines show model spectra for the Be star and donor; black line shows their sum after accounting for extinction (Ay = 0.8 mag). These are not a fit 
to the full observed SED, but are scaled to match the blue optical spectrum (Figure 4). Bottom panel shows the ratio of the observed fluxes to the predicted 
bandpass-integrated magnitudes. There is definite evidence of IR excess from the disk, which contributes ~ 20% of the total light in the H band, and more at 


longer wavelengths. 


affected by CNO processing. We therefore added the helium, carbon, 
nitrogen, and oxygen abundances as additional labels to be fit. 

We implicitly assume that the two stars have the same surface 
abundances. In practice, the metallicity and CNO abundances are 
well-constrained only in the donor, since metal lines are rotationally 
smeared out in the Be star, and the He abundance is constrained 
only in the Be star, due to the donor's cooler temperature. We mask 
the Balmer and Brackett lines during fitting, since they are contami- 
nated by emission from the disk. We also mask the two helium lines 
with obvious emission, and regions that are affected by tellurics or 
interstellar absorption (Figure 5). 

Among the largest sources of uncertainty in fitting the spectrum 
is the micrturbulence of the donor, which is covariant with the flux 
ratio (and thus with Rdonor and Teff, donor). This is illustrated in Fig- 
ure 10, which shows best-fit models for the PEPSI spectrum for three 
different values of Viy;p. Increasing Vturp increases the equivalent 
width of most metal lines. This has a (qualitatively) similar effect to 
increasing the flux ratio or decreasing the donor's effective tempera- 
ture. The best-fit microturbulence is vy = 2.5 km sl, larger than 
typical value of ~ 1 km s^! for main-sequence stars. 

The spectra have per-pixel SNR between 500 and 2000, so full 
spectral fitting produces unreasonably small formal uncertainties 
(e.g., a few Kelvin in Ta). These dramatically underestimate the 
true uncertainties — to the point of not being particularly useful — be- 
cause they do not account for systematics in the model spectra (e.g., 
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imperfect linelists, assumption of LTE, 1D treatment of convection), 
systematics in the spectral fitting (e.g., assuming most metal abun- 
dances track Fe, imperfect normalization), or systematics in the data 
(e.g. imperfect masking of tellurics and DIBs, bad pixels, unrecog- 
nized contamination from emission lines, etc). This is a perennial 
issue for full spectral fitting in the high-SNR regime (e.g. Ness et al. 
2015). To obtain a more conservative estimate of the uncertainties, 
we varied the microturbulent velocity by +1 km s^! around its fidu- 
cial value of 2.5 km s^! (i.e., fixing itto 1.5 and 3.5 km s7!) and re-fit 
the other parameters. The + uncertainties we report are half of the 
minimum-to-maximum range across the three fixed values of vq, 
with the reported best-fit values corresponding to vturb = 2.5 kms !. 
This approach leads to serviceable uncertainties that are more be- 
lievable than the formal fitting errors, while still capturing the most 
important sources of uncertainty and the covariances between pa- 
rameters. 


4.2 Be star temperature, radius, and mass 


Figure 4 shows the flux-calibrated spectrum and the predicted spec- 
tra of the best-fit model. The Be star's effective temperature and 
radius are only weakly sensitive to the properties of the donor, which 
contributes € 20% of the blue optical light. 

We compare the inferred parameters of the Be star to MIST evolu- 
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Figure 10. PEPSI spectrum and best-fit Kurucz model spectra for three different values of v5, the microturbulent velocity. In all panels, we fix the parameters 


of both components to the values listed in the top panel. We fix vq to a different value in each panel and solve for the Rdonor value that best fits the observations 
given this vturb. Increasing Rgonor increases the fraction of the flux contributed by the donor. Conversely, increasing vturb makes most of the donor's lines deeper. 
There is thus a significant covariance between Vturb and Rdonor. For the best-fit vy = 2.5 km s^! and Rgonor = 5.8 Ro, the donor contributes 21% of the total 


flux at 5200 Á. 


tionary models (Choi et al. 2016) for solar-metallicity main-sequence 
stars in Figure 11. To estimate the mass of the Be star, we construct a 
grid of evolutionary tracks with a mass spacing of 0.05 Mo, drawing 
500 samples spaced uniformly in time from each track. We then re- 
tain samples with probability proportional to the likelihood of their 
effective temperature and radius given the measured constraints on 
the Be star's atmospheric parameters, approximating the likelihood 
function as a 2D Gaussian. This effectively selects points from the 
evolutionary tracks that are close to the Be star in temperature and 
radius. We use the resulting samples to estimate the mass and lu- 
minosity of the Be star, and to propagate forward the uncertainty in 
other parameters that depend on the Be star's properties. This yields 
a Be star mass Mpe = 5.3 + 0.6Mo. Implicit in this modeling is 
the assumption that the Be star falls on a mass/radius/temperature 
relation similar to normal main-sequence stars (see Section 4.4). 


We note that we could also calculate Mp. directly from the mea- 


sured surface gravity and radius. This leads to Mge = 7.7 € 4.3 Mo. 
This value is consistent with our inference from the evolutionary 
models, but its uncertainty is large due to the uncertainty in log g. 


4.3 Donor temperature, radius, and mass 


The donor effective temperature and radius are also well-constrained 
by the spectra. The radius constraint comes from the normalization 
of the flux-calibrated spectrum and constraint on the optical flux ratio 
from the high-resolution spectrum. We find Rgonor = 5.82 + 0.45 Ro 
and Teff, donor = 6900 + 250K, with the flux ratio the domi- 
nant source of uncertainty. The spectroscopic surface gravity is 
log(gdono;/cm s?) = 2.8 + 0.15. 

The mass of the donor can be constrained from the fact that it 
must fit within an orbit with Porp = 5.47 days; i.e., the donor must 
be massive enough that its observed radius is not larger than its 
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Figure 11. Measured parameters of the Be star in HD 15124 compared to 
MIST evolutionary tracks for single stars. Effective temperature and surface 
gravity are constrained spectroscopically; radius is constrained by the spectral 
energy distribution. The observed parameters imply a mass Mge = 5.3 + 
0.6Mo. 


Roche lobe radius. Taking the expression for the Roche lobe radius 
from Eggleton (1983) leads to a lower limit on Mgonor for any given 
Raonor that depends only very weakly on Mge. This constraint is 
shown in Figure 12, and is essentially a reflection of the fact that all 
Roche lobe-filling stars at a given Pot have nearly the same density, 
irrespective of their companion mass. 

With a lower limit on the donor mass in place from the Roche 
lobe constraint, an upper limit can be set in multiple ways. First, 
the spectroscopic log g = 2.8 + 0.15 and radius constraints imply 
Maonor = (0.77 + 0.32) Mo, or Maonor < 1.09 Mo, not much larger 
than the Roche lobe lower limit. Second, the amplitude of ellip- 
soidal variation depends strongly on the Roche lobe filling factor 
(e.g. El-Badry et al. 2021d, their Figure 10). Based on the sec- 
ondary minima present in the TESS light curve (Figure 6), we find 
Rdonor/RRochelobe 2 0.93. This translates to a donor mass only 
~ 25% larger than when Rgonor/RRochelobe = 1. For Rdonor = 5.82, 
requiring that 0.93 < Rgonor/RRochelobe € | constrains the donor 
mass to 0.92 < Mdonor/ Mo < 1.15. 

Since the light curve and spectroscopic log g suggest that the donor 
is close to Roche lobe-filling, we suspect that mass transfer is still 
ongoing and the donor is semi-detached. This naturally explains 
the fact that the Be star has a circumstellar disk while rotating at 
a lower fraction of its critical rotation velocity than most Be stars 
(Section 4.5): it is still being spun-up by accretion. Moreover, all the 
evolutionary models we construct to match the system (Section 5) 
are still mass-transferring when they have temperature and surface 
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Figure 12. Donor mass, donor surface gravity, and binary semi-major axis 
predicted for different donor radii in HD 15124 if the donor fills its Roche 
lobe, as we argue it must. Shaded region shows the range of plausible donor 
radii, given our fits to the SED and spectrum. Hatched region shows the 
observational log g constraint. The donor mass and surface gravity are nearly 
independent of the mass of the Be star for any fixed donor radius; these 
quantities depend only on the (known) orbital period and the donor mass 
and radius. If the donor did not fill its Roche lobe, its mass, surface gravity, 
and semi-major axis would all be /arger than these predictions. If it does, 
our constraint on Rdonor implies Maonor = 0.92 + 0.22 Mo, log Zdonor = 
2.87 + 0.04, and a = 24+ 0.8 Ro. 


gravity similar to the donor in HD 15124. We thus proceed under the 
assumption that the donor exactly fills its Roche lobe. This leads to a 
donor mass constraint of Mgonor = 0.92 + 0.22 Mo. The uncertainty 
is dominated by the uncertainty on Rgonor, and thus would be only a 
factor of ~1.5 larger if we did not treat the system as semi-detached. 


4.4 Inclination 


From the measured RVs of the donor, we calculate a mass function 
f(Mpe) = 0.23 + 0.01 Mo. This quantity represents the absolute 
minimum mass of whatever object the donor is orbiting. In particular, 
itis the mass that that object would have if the donor were a massless 
test particle and the system were viewed edge-on. For a finite donor 
mass or lower inclination, the implied mass is higher. 

We proceed under the assumption that the donor is orbiting the Be 


star, and interpret f (Mpe) as a constraint on the Be star's mass. This 
assumption is not guaranteed to hold, since we have not measured 
RVs for the Be star to demonstrate the two stars are orbiting each 
other. It is in principle possible that the donor is orbiting another faint 
object (say, a white dwarf) and the Be star is a distant tertiary or a 
chance projection. However, we consider both scenarios unlikely. The 
chance projection scenario can be rejected because the probability 
of a close alignment of two bright stars with less than one arcsec 
separation is very small: following the approach used by El-Badry 
et al. (2021c), we find a chance alignment probability of x 1076. The 
triple scenario with a faint close companionis less unlikely, but would 
not explain the observed reflection effect in the light curve (Figure 6). 
We thus assume the donor is orbiting the Be star. We discuss the triple 
scenario with an outer Be star further in Appendix C. 

The observed mass function constrains the binary to a 2- 
dimensional surface in the space of donor mass, companion mass, and 
orbital inclination. This is illustrated in the bottom panels of Figure 8, 
which show the implied Be star mass for various combinations of 
donor mass and inclination. The relatively low mass function implies 
either a low companion mass or a low inclination (close to face-on). 
Given our constraints on the masses of both stars, the mass function 
implies an inclination 7 = 23.2 + 1.5 degrees. We also consider the 
possibility that our spectroscopic constraint on the companion mass 
could be catastrophically wrong (e.g., if the Be star were a distant 
tertiary to an inner binary). In this case, the mass function would 
imply a companion mass of 0.92 + 0.02Mo for an inclination of 90 
degrees, or 1.15 x 0.02Me for an inclination of 60 degrees. However, 
we consider such a scenario unlikely (Appendix C). 

Only about 8% of randomly oriented orbits have inclinations as 
low as 23 degrees. It is thus worth considering whether there could be 
any serious systematic error affecting the inclination measurement. 
The inferred inclination depends mainly on the mass function and 
the mass of the Be star. The mass function is quite secure. The Be 
star's mass was inferred by comparing its spectroscopic temperature 
and radius to evolutionary models (Figure 11). This inference treats 
the Be star as a main sequence star and is thus somewhat less secure, 
because the Be star could in principle be out of equilibrium due to 
recent accretion. However, our MESA models (Section 5) suggest 
that departures from equilibrium for the accretor are modest, such 
that the accretor temperature and radius remain similar to those of 
a main-sequence star of the same mass. Even if our estimate of the 
Be star's mass were too large by a factor of 2 (which seems very 
unlikely), the implied inclination would remain low, i ~ 33 deg. We 
conclude that the inclination is genuinely low. 


4.5 Rotation velocities 


We find projected rotation velocities v sini = 20 + 3kms~! for the 
donor and v sini = 95 + 10 km s^! for the Be star. The constraint for 
the donor comes primarily from metal lines in the PEPSI spectrum, 
while that for the Be star is from helium lines. If we assume the 
spin axes of both stars are aligned with the binary’s orbital angu- 
lar momentum vector — a reasonable assumption for a close mass- 
transferring system — then the constraint on the orbital inclination 
allows us to translate these into a synchronization parameter for the 
donor, vrot/(27R/Porp) = 0.94 + 0.17; i.e., the donor is probably 
tidally synchronized. This validates our inferred inclination: a near 
Roche-filling star is expected to be tidally synchronized, and the mea- 
sured v sini would not correspond to synchronous rotation if it were 
not for the low inclination. 

The implied critical rotation fraction for the Be star is vrot/Verit = 
0.60+0.08. This is a factor of 2-3 larger than typical values for normal 
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Table 2. Photospheric abundances. We assume both stars have the same 
surface abundances. In practice, the constraint on He comes primarily from 
the Be star, while the constraints on all other elements are from the donor. 
Other metals track Fe, assuming the solar abundance pattern from Asplund 
et al. (2009). Uncertainties in the 4th column also apply to the 1st column. 


Element log(n/ntot) —_ log(n/niot)o log (er 


(n/niot) o 
He -0.51 -1.11 0.6 x 0.1 
C -5.91 -3.61 -2.3 40.2 
N -3.41 -4.21 0.8 + 0.2 
[9) -4.15 -3.35 —0.8 + 0.2 
Fe -4.49 -4.54 0.05 x 0.1 


(non-emission line) B stars (Abt et al. 2002; Huang et al. 2010), but 
is also lower than the typical values of « 0.8 — 0.9 for classical Be 
stars (Townsend et al. 2004). This suggests that the star has been spun 
up by accretion, but has not yet reached critical rotation. 


4.6 Be star disk 


The velocity separation between the two peaks of the emission lines 
ranges from 250 to 350 km s^, while the separation between the 
line wings ranges 500 to 600 km s7! (Figure 5). If we assume that 
the emission is optically thin, then the separation between the peaks 
traces the Keplerian velocity at the outer edge of the disk, since 
for typical emissivity profiles, that is where most of the emission 
originates (e.g. Smak 1969). Given the inferred inclination, the peak- 
to-peak velocity separation would correspond to orbital disk veloc- 
ities of Vouter % 380 70km s7!, corresponding to an outer disk 
radius of Router = GMpe/V> ster = (7.0 + 3.5) Ro, or in terms of 
the orbital separation, Router = (0.29 + 0.14) a. This is consistent 
with the maximum stable disk radius of Rmax ~ 0.42a for a binary 
with mass ratio « 0.2 (Paczynski 1977), so the emission may trace 
material out to the largest stable streamline. We caution, however, 
that the emission is not necessarily completely optically thin, and the 
above does not hold if the line profile is significantly broadened by 
non-coherent scattering (i.e., absorption and re-emission at a differ- 
ent wavelength). Indeed, scattering must contribute to broadening at 
least the line wings, since their velocity separation would correspond 
to (2.1 + 0.6Ro) - inside the Be star — if it were purely Keplerian. 


4.7 Chemical abundances 


The surface abundances (Table 2) reveal significant enhancement of 
helium and nitrogen, and depletion of carbon and oxygen. Figure 13 
compares cutouts of the high-resolution spectra to the best-fit model 
(cyan) and a model with the same atmospheric parameters and flux 
ratios but Solar abundances. The helium abundance is constrained by 
several strong lines in the PEPSI spectrum, and the carbon abundance 
by many lines in both the PEPSI and APOGEE spectra. The oxygen 
and nitrogen abundances are determined by a smaller number of 
weaker lines but still appear well-constrained. 

The best-fit model has a helium abundance by number that is a 
factor of 4 (0.6 dex) larger than the solar value. This corresponds to a 
surface helium mass fraction Yye = 0.63 + 0.08. Carbon is depleted 
by 2.3 dex (more than a factor of 100), while oxygen is depleted by 0.8 
dex and nitrogen is enriched by 0.8 dex. The carbon-to-nitrogen ratio 
is thus a factor of ~ 1000 larger than the solar value. The abundances 
of other elements appear normal. 

These abundances manifestly imply that the material currently 
on the surface of both stars was previously processed in the CNO 
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cycle, which converts hydrogen to helium using C, N, and O as 
catalysis. In stars with masses of a few to a few tens of solar masses, 
CNO burning leads to a chemical equilibrium of C, N, and O in the 
convective core while the star is on the main sequence (e.g. Maeder 
et al. 2014). Once the outer envelope of the donor is stripped off, this 
material is exposed in the photosphere. HD 15124 thus provides a 
rare opportunity to effectively look inside the core of an intermediate 
mass star. As we show in Section 5, the current surface abundances 
are in good agreement with theoretically predicted values in such a 
scenario. 


4.8 Galactic orbit 


To investigate the past trajectory of HD 15124, we used its Gaia 
astrometry and inferred center-of-mass radial velocity as starting 
points to compute its orbit backward in time for 500 Myr (though we 
note that the evolutionary age of the system is only 100-200 Myr) 
using galpy (Bovy 2015), adopting the Milky Way potential from 
McMillan (2017). The result is shown in Figure 14. The orbit appears 
typical of a thin-disk star, with excursions above the disk midplane 
limited to +50 pc. We verified that this conclusion is robust to changes 
in the assumed potential and to varying the astrometry and RV within 
the observational uncertainties. 


5 EVOLUTIONARY HISTORY OF THE HD 15124 SYSTEM 
5.1 BPASS library search 


To understand the possible formation pathways and future evolution 
of HD 15124, we began by searching the BPASS (v2.2; Eldridge 
et al. 2017; Stanway & Eldridge 2018) library of binary evolution 
calculations for models that go through a phase similar to the binary's 
current state. We searched the Z = 0.02 grid for models that at any 
point in their evolution satisfy the following constraints: 


5900 < Teff, aono; / K < 7900 
4 < Rgonor/ Ro < 8 

4 « M3/ Mo < 10 

0 < Porb/days < 10 


Here Mp is the mass of the initially less-massive star. This search 
yielded 18 models, with initial primary masses ranging from 3.5 to 
5 Mo, initial mass ratios ranging from 0.5 to 0.9, and initial periods 
of 1.6 to 2.5 days. Their evolutionary history is qualitatively quite 
similar to the models explored to explain the HR 6819, LB-1, and 
NGC 1805 BH1 systems (e.g. Eldridge et al. 2020; Bodensteiner 
et al. 2020b; El-Badry & Quataert 2021; El-Badry & Burdge 2021; 
Stevance et al. 2021). In all the selected models, mass transfer begins 
while the donor is still on the main sequence (“case A’). While these 
models do come close to the observed properties of HD 15124, they 
all predict donor masses and radii on the high side of the observed 
values (median Mgonor ¥ 1.7 Mo and Raonor & 7.1 Ro). The surface 
abundances of the donors in these models suggest that they are less 
stripped than HD 15124 (median Yye = 0.31, only slightly enriched 
compared to the initial value of 0.28). 

Models with lower donor masses and higher Ype do exist in the 
BPASS library, but these have longer periods. This is a consequence 
of the way mass and angular momentum transfer are implemented 
in the BPASS models: mass transfer is assumed to be fully conser- 
vative, unless it occurs on a timescale shorter than the accretor’s 
thermal timescale (e.g., M > Maccretor/ tKH, accretor- Where fkg is 


MNRAS 000, 1—30 (2022) 


the Kelvin-Helmholtz time). This leads to rapid widening of the or- 
bit once the donor is the lower-mass star (e.g. Equation 1 below). 
To investigate the evolution of the system in more detail and ex- 
plore models with different treatments of angular momentum loss, 
we calculated additional models. 


5.2 MESA calculations 


We calculated a small grid of binary evolution models using MESA 
(Modules for Experiments in Stellar Astrophysics, version r15140; 
Paxton et al. 2011, 2013, 2015, 2018, 2019). We primarily explored 
initial conditions similar to those identified by our BPASS search, 
but also experimented with varying the initial conditions and input 
physics. Unlike the models in the BPASS library, which only follow 
the primary star with detailed calculations, MESA simultaneously 
solves the 1D stellar structure equations for both stars, while ac- 
counting for mass and angular momentum transfer using simplified 
prescriptions. 

We assumed an initial composition of X = 0.7, Y = 0.28, Z = 0.02 
for both stars. Opacities are taken from OPAL (Iglesias & Rogers 
1996) at log(T/K) > 3.8 and from Ferguson et al. (2005) at 
log(T/K) « 3.8. Nuclear reaction rates are from Angulo et al. (1999) 
and Caughlan & Fowler (1988). We use the pp. and. cno. extras 
nuclear network, which includes 12H, 3.4He, "Li, 7Be, 8B, 14n, 
14,160, 19F, 18,19,20NŅe, and 22,24Mg. We use the photosphere 
atmosphere table, which uses atmosphere models from Hauschildt 
et al. (1999) and Castelli & Kurucz (2003). 

The MESAbinary module is described in Paxton et al. (2015). We 
used the evolve_both_stars inlists in the MESA test suite as a 
starting point for our calculations, and most inlist parameters are set 
to their default values. Both stars are evolved simultaneously, with the 
response of the secondary to accretion taken into account. Roche lobe 
radii are computed using the fit of Eggleton (1983). Mass transfer 
rates in Roche lobe overflowing systems are determined following the 
prescription of Kolb & Ritter (1990). The orbital separation evolves 
such that the total angular momentum is conserved when mass is lost 
through winds or transferred to a companion, as described in Paxton 
et al. (2015). 

We experimented with both conservative and non-conservative 
mass transfer. As with the models from the BPASS library, we were 
unable to match the current properties of HD 15124 (orbital period, 
donor mass and radius, and significant surface helium enrichment) 
with models in which mass transfer is fully conservative (meaning 
that all mass lost by the primary is accreted by the secondary). 
The “problem” is the short orbital period. When mass transfer is 
conservative, there is a simple relation between the initial and current 
orbital periods and mass ratios (e.g. Soberman et al. 1997): 


3 
M donor, init M accretor, init 
Maonor ( t ) Maccretor ( t ) 


Here Porb(t), Maonor(t), and Maccretor(t) are the current orbital pe- 
riod and donor and accretor masses, while values with the subscript 
"jnit" are their values at the onset of mass transfer. Given the cur- 
rent Mdonor & 0.92 Mo and Maccretor & 5.3 Mo, the total initial 
mass would need to be about 6.2M if mass transfer were fully 
conservative. To avoid the formation of a common envelope or con- 
tact binary, the initial masses are required to have been relatively 
similar; e.g., Maccretor, init/ Mdonor,init 2 0.6. Taking, as a repre- 
sentative example, Porb, init = 2 days and Maonor, init = 3.5Mo and 
Maccretor, init = 2-7 Mo, Equation 1 yields a current period Porp ~ 15 
days, significantly longer than the observed value of 5.5 days. In 
order to match the observed period with these masses, Equation 1 


Pow (t) z Porb, init x (1) 
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Figure 13. Spectral cutouts highlighting evidence of CNO processing. Top two rows show PEPSI spectra; bottom row shows APOGEE spectra. In all panels, 
black line shows the observed spectra, and colored lines show composite model spectra with the parameters listed in Table 1. Helium lines are primarily from 
the Be star; all narrow lines are primarily from the donor, with continuum dilution from the Be star. Red line shows the best-fit model, the abundances of which 


are listed in the upper left panel. Cyan shows solar abundances, and magenta (bottom panels) shows other choices of N and C abundances. The spectra reveal 
strong enhancement of helium and nitrogen at the expense of carbon and oxygen, a sign of CNO processing. This suggests that the material currently on the 


surface of both stars was previously inside the convective core of the donor. 


requires a shorter initial period. But in this case, the donor would 
overflow its Roche lobe early in its main-sequence evolution, leading 
to qualitatively different evolution (Section 5.3). 


When mass-transfer is not conservative, the orbital period evolu- 
tion depends on the angular momentum of the material that escapes. 
MESABinary implements the “afy6” formalism (see Huang 1963; 
van den Heuvel 1994; Soberman et al. 1997; Tauris & van den Heuvel 
2006) for non-conservative mass loss, in which a fraction ayy of the 
mass lost by the donor is ejected in a fast wind from the donor, a 
fraction By is ejected in a fast wind from the accretor, and a fraction 


Omi is lost from a circumbinary ring with radius Nr X a. In our 
fiducial model, we fix ayy, = Bi = 0, Omi = 0.2, and YML = V2; 
that is, we assume 20% of the mass leaving the donor is eventually 
lost from the binary through a circumbinary ring of radius 2a, while 
the rest is transferred to the companion. Compared to the case of con- 
servative mass loss, this tends to reduce the binary’s orbital period, 
because the mass that is lost from the circumbinary ring has higher 
specific angular momentum than the binary itself. 


In reality, «ML, BML; YML, and Opp are expected to be functions 
of the mass transfer rate and other instantaneous properties of the 
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Figure 14. Galactic orbit of HD 15124, computed backwards for 500 Myr 
McMillan (2017). Red and yellow stars show the present-day location of HD 


binary, and thus should vary with time. Ab-initio constraints on these 
parameters are scarce, and we do not attempt to infer values for them 
based on one system. For our purposes, the choice of ôM, = 0.2 
and yML = v2 should be regarded as an effective model of a more 
complicated mass transfer process; likely, one of several possible 
combinations of «ML, BML; YML, and Om that can produce the right 
amount of angular momentum loss to match the observed period. We 
discuss how varying the parameters of the MESA model from our 
fiducial values changes the calculation in Section 5.3. 


Figure 15 shows the evolution of the initially more massive star 
in a MESA model that matches the present-day properties of HD 
15124, which we adopt as our fiducial model. The initial component 
masses are 3.75 and 3.0 Mo, and the initial period is 2.4 days. Mass 
transfer begins as the primary is crossing the Hertzsprung gap and 
initially proceeds on its thermal timescale, leading to a maximum M 
of x 107? Mo yr! . Here the donor rapidly cools, and its luminosity 
falls by more than a factor of 10. During this phase, the deeper layers 
of the donor are expanding to adjust to the decrease in total mass. The 
energy needed for this expansion comes at the expense of the surface 
luminosity, such that the nuclear luminosity significantly exceeds the 
surface luminosity. 


Once the donor has lost most of its envelope, the mass transfer 
rate falls to M = 1077 Mo yr! and the donor begins to expand. By 
now the donor has become the lower-mass star, so mass loss widens 
the orbit. This leads to a semi-detached system whose evolution is 
governed by expansion of the evolving donor. During this period, 
the model resembles HD 15124, with Mgonor 7 0.8 Mo, Tem donor ~% 
7000 K, and log(gdonor) « 2.8. 


The evolution of the model during this period is shown in more 
detail in Figure 16. The period during which our APOGEE search 
would be sensitive to the donor lasts « 1 Myr, after which the donor 
becomes too hot to be detectable via its metal lines in the IR (Ap- 
pendix B). The binary remains semi-detached for another ~ 1 Myr, 
until the donor’s hydrogen envelope mass is reduced to ~ 0.1 Mo. At 
this point, the donor begins to contract and heat up, eventually set- 
tling at a radius of z 0.2 Ro and effective temperature of ~ 35, 000 K, 
typical values for a sdO/B star (Heber 2016). 

The luminosity of the model during the observed phase and the first 
= 2 Myr of contraction is powered by shell burning of the remaining 
hydrogen envelope. Core helium burning begins just as the donor 
begins to contract and dominates the luminosity after ~ 2 Myr. The 
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from the measured astrometry and center-of-mass RV using the potential from 
15124 and the Sun. The orbit is typical for a young star in the Galactic thin disk. 


lifetime of the subsequent core helium burning phase would be of 
order 100 Myr, but in most of the models we explore, the companion 
finishes its main sequence evolution and overflows its Roche lobe 
first. This is likely to lead to a period of common envelope evolution, 
but we stop the calculation here. 

In Figure 15, we also compare the model to the donor stars in 
LB-1, HR 6819, and NGC 1850 BHI. These systems can all be de- 
scribed by a qualitatively similar evolutionary scenario in which an 
intermediate-mass primary overflows its Roche lobe near the end of 
its main-sequence evolution, goes through an Algol phase, and even- 
tually contracts to become an sdO/B star. The details are, however, 
different in each system. LB-1 and HR-6819 have significantly longer 
orbital periods (79 and 40 days) and are thus better described by mod- 
els with more conservative mass transfer. HD 15124 has a shorter 
orbital period than any known classical Be star binaries. Several 
Algol-type binaries with qualitatively similar evolutionary histories 
and similarly short periods are known (e.g. Tupa et al. 2013; Ros- 
ales et al. 2021), but the donors in these systems likely have higher 
masses and are somewhat less evolved than HD 15124. The orbit 
of HD 15124 is expected to widen as mass transfer continues. The 
fiducial MESA model predicts a final period of 12 days, but the true 
value may be larger because mass transfer is expect to become more 
conservative as M decreases, an effect not included in the model. 

Figure 17 shows the evolution of surface abundances of the donor 
in the MESA model as a function of orbital period. These are ini- 
tially solar but reveal increasing levels of CNO processing as the 
period increases and deeper layers of the donor's interior are ex- 
posed. Figure 18 shows the interior abundance profiles of the donor 
at several snapshots of its evolution. CNO ratios at the surface remain 
at the solar value until initial Roche lobe overflow, when the core is 
already heavily processed. Helium and oxygen are respectively en- 
riched and depleted only in the central x 1 Mo, while nitrogen and 
carbon are enriched and depleted out to 2 Mo. The observed surface 
abundances thus provide an additional constraint on how much of 
the donor's outer envelope has been stripped. 


5.3 Varying model parameters 


Here we discuss how the behavior of the MESA models changes as 
we vary the initial conditions and model ingredients. Some represen- 
tative models are shown in Figure 19. 
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Figure 15. MESA binary evolution calculations for a system similar to HD 15124. Let panels show the evolution of the primary (i.e., the donor) in a binary 
with initial masses of 3.75 and 3 Mo and initial period of 2.4 days. During its main-sequence lifetime, the primary evolves from the zero-age main sequence 
(ZAMS) to the Hertzsprung gap (HG). Roche lobe overflow (RLOF) occurs as the star is crossing the HG and is followed by a brief (~ 3 x 10° years) period 
of thermal-timescale mass transfer. The star cools at near-constant radius during this period, leading to a factor-of-40 drop in luminosity. The donor reaches 
Toe ~ 7000 K and log g ~ 2.8 (similar to HD 15124) at Py = 6 — 7 days, when its mass loss rate is M ~ 1077 Mo yr~!. Mass transfer ceases after ~2 Myr, 
when the donor has lost most of its envelope. The donor then contracts and heats up, settling as a core helium burning sdO/B star. During this contraction, the 
model’s parameters are similar to those of the bloated stripped stars in HR 6919, LB-1, and NGC 1805 BH1, systems with similar evolutionary histories. The 
calculation ends after another ^50 Myr, when the Be star evolves and overflows its Roche lobe. 
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Figure 16. Time-evolution of the fiducial MESA model, zoomed-in on the 
period of interaction. Mass-transfer begins at t = 0 and ends at t = 2.5 Myr, 
after which the donor heats up and contracts. Most of the mass loss occurs 
in the first 0.5 Myr. Shaded region indicates when the donor temperature 
and surface gravity are similar to HD 15124, such that our APOGEE search 
would detect the companion. Third panel shows the donor luminosity (black) 
and luminosity of hydrogen burning (red), helium burning (cyan) and grav- 
itational contraction (magenta). The luminosity is currently dominated by 
shell hydrogen burning but is predicted to become dominated by core helium 
burning within a few Myr. Bottom panel shows the fraction of the total flux 
contributed by the donor at optical, UV, and NIR wavelengths. 
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Figure 17. Surface abundances of HD 15124 and the MESA model from 
Figure 15. All panels show mass fractions, and dashed horizontal lines show 
Solar values. The observed abundances of helium, carbon, nitrogen, and 
oxygen are all highly discrepant from the solar abundance pattern: helium 
and nitrogen are strong enriched, and carbon and oxygen are depleted. The 
model reproduces the observed abundances reasonably well. In the model, 
the current surface of the donor was previously inside the donor’s convective 
core during CNO burning. 


(i) Larger or smaller initial period: A moderately longer initial 
period (3 € Pinit/day < 20, with all else fixed) causes mass transfer 
to begin farther up the subgiant branch, when the donor is larger. The 
model’s evolution is qualitatively similar to in our fiducial model, 
but the current orbital period is longer than observed. At even longer 
initial periods (20 € Pinit/day x 150), mass transfer begins when 
the donor is a giant, presumably leading to an episode of common 
envelope evolution (which cannot be properly followed in our cal- 
culations). This could potentially produce a short period as in HD 
15124 but likely would not end with a donor that is still Roche lobe- 
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Figure 18. Profiles of CNO mass fractions in the donor, for the MESA model shown in Figure 16. The four panels show the zero-age main sequence, onset 
of Roche lobe overflow, 1 Myr later (corresponding roughly to the current state of HD 15124), and core helium burning. The horizontal axis changes between 
panels as mass is removed from the donor. Vertical line in panel 2 marks the approximate current surface of the donor. The surface of the star today fell within 
the convective core of the donor at RLOF, whose composition is heavily affected by CNO processing. 


filling. Even longer initial periods can produce stable mass transfer 
when the donor is an AGB star (stabilized by its larger core mass), 
but these also produce periods much longer than observed. 

A shorter initial period (Pig < 2 days) causes mass transfer to 
begin while the donor is still on the main sequence. Mass transfer is 
initially similar to in the fiducial model. However, because the donor 
has not yet terminated its nuclear evolution, it shrinks back within 
its Roche lobe once it becomes the lower-mass star, and further 
evolution occurs on a (now longer) nuclear timescale. Eventually, 
the now-more-massive secondary evolves off the main sequence and 
overflows its Roche lobe. Most of the BPASS models we explored 
follow this type of evolution, but they do not track the response of 
the secondary to accretion and thus do not capture its more rapid 
subsequent evolution. 

(i) Higher or lower primary mass: A higher initial primary mass 
is possible if (a) the secondary mass was also higher (because a 
too-unequal mass ratio leads to formation of a common envelope; 
see below) and (b) the mass transfer efficiency was lower, to avoid 
producing a too-massive Be star. Under these conditions, a higher 
initial primary mass generally leads to a higher current mass for the 
bloated core. A significantly lower initial primary mass is ruled out, 
since the total mass of the binary today must exceed ~ 6 Mo. 

All the MESA models we found plausibly compatible with the 
current properties of HD 15124 produce helium cores with masses 
0.5 € Mgecore/ Mo X 1, such that the donors will become core 
helium burning sdO/B stars. Lower-mass analogs of the system pro- 
duce donors that will contract directly to become helium white dwarfs 
(e.g. El-Badry et al. 2021b; Miller et al. 2021). Higher-mass systems 
produce more massive and hotter sdO/B stars (Gótberg et al. 2017), 
which eventually can explode and leave behind a Be X-ray binary 
(e.g. Townsend & Charles 2020; Klement et al. 2021a). 

(iii) Larger or smaller mass ratio: Mass transfer is less stable 
(i.e., more likely to result in a common envelope) when the ratio 
Maonor/Maccretor is large. Models with highly unequal initial masses 
(M2/M, < 0.6) can in principle be evolved with MESA to produce 
systems like HD 15124 and similar binaries at more evolved stages 
(e.g. Bodensteiner et al. 2020b), but this entails evolving the model 
through a contact phase, which is not reliably captured by the models. 
Mass ratios closer to 1 can produce qualitatively similar evolution 


to our fiducial model. For initial mass ratios close to 1, the duration 
of the predicted Be + sdO/B phase becomes shorter, because the 
secondary leaves the main sequence sooner. 

(iv) Changes in mass transfer prescription: If mass loss is con- 
servative, the orbital period increases too rapidly to reproduce the 
observed value (Equation 1). This can produce wider systems like 
LB-1 and HR 6819. If mass loss is less conservative or removes more 
specific angular momentum than in our fiducial model, the orbital 
period will be shorter, and the donor will be fully stripped before 
it can ascend up the subgiant branch. This can produce a compact 
mass-transferring system like NGC 1850 BHI. 


6 IMPLICATIONS FOR BE STAR FORMATION 
CHANNELS 


Once the donor in HD 15124 becomes an sdO/B star, it will be 
outshone by the Be star and very difficult to detect at any wavelength 
(bottom panel of Figure 16). Its mass at that time is predicted to be 
z 10 times lower than that of the Be star, so the expected RV shifts of 
the Be star will be only a few km s^! — too low to reliably detect in a 
Bestar with high rotation velocity and disk-driven spectral variability. 
The lifetime of the Be + sdO/B phase is expected to be much longer 
than that of the current phase when the donor is detectable. This 
suggests that there are very likely many Be + sdO/B binaries within 
the sample of Be stars from which we selected HD 15124 that have 
thus-far undiscovered sdO/B companions. 

We detected Naetecieg = 1 bloated stripped companions among 
Nge = 297 Be stars. From this, we can infer fgetectable» the fraction 
of Be stars with bloated stripped companions detectable with our 
search. From Bayes' theorem, 


P Cfüetectablel Naetected) © p (Naetected| fdetectable) P (fdetectable) - 
(2) 


Here p (Naetected| fdetectable) is the likelihood of detecting Naetected 
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Figure 19. Effects of varying MESA model parameters. Top panel shows the fiducial model. Second panel shows a longer initial period, which causes mass 
transfer to begin later, such that significant helium enrichment is not predicted until longer periods than observed. Third panel shows a shorter initial period. 
In this case, mass transfer begins before core hydrogen exhaustion in the donor but then stalls until the now-more-massive secondary evolves and overflows its 
Roche lobe. Third panel shows a larger value of yy, such that the specific angular momentum of matter lost from the binary is higher. This leads to more 
rapid stripping and a shorter orbital period. Bottom panel shows fully conservative mass transfer. In this case, the donor does not reach the observed degree of 
stripping as quantified by Ype until a much longer orbital period. 


MNRAS 000, 1—30 (2022) 


stripped companions given faetectable- For a Poisson processes, this 
is given by 


N, 
Cfaetectabie NBe) ads) g^ actectable NBe (3) 


P (Naetected| faetectable) = N. 1 
detected: 


For the prior, p ( faetectapie). We assume a uniform distribution be- 
tween 0 and 1. The proportionality factor in Equation 2 is then set by 
the normalization condition. 

To relate fdetectable tO fstripped» the total fraction of Be stars that had 
or have stripped companions, we need an estimate of Tge, the lifetime 
of the Be phase, and Tgetectable, the lifetime of the phase during which 
a typical companion would have been detectable with APOGEE. For 
simplicity, we assume Tgpe = 50 Myr and Tgetectapile = 1 Myr, as 
suggested by our evolutionary models for HD 15124, so that 


stripped E TBe — 50. (4) 


detectable — "detectable 


Figure 20 shows the resulting posterior distribution of fstripped- 
The distribution is broad — due entirely to Poisson uncertainties 
— but suggest that a significant fraction of Be stars have stripped 
companions. The median and middle 68% of the distribution is 
stripped = 0.28+9:27; i.e., discovery of one stripped companion im- 
plies that (10-60) % of Be stars have previously had a bloated stripped 
companion. In most cases, these companions should survive as faint 
white dwarf, sdO/B, and neutron star companions. 

We note that HD 15124 is in the brightest 20% of the Be stars in 
initial sample (Figure 1), raising the question whether there could 
be additional similar stripped companions in the sample that were 
missed. If we had detected a 2nd similar object, the above constraint 
would change to fstripped = gd. However, the stripped com- 
panion in HD 15124 is not near the limits of our search's sensitivity 
(Appendix B). 

This estimate relied on numerous approximations. Perhaps most 
significantly, our presumed lifetimes TBe and Tdetectable are based on 
models for HD 15124 and will vary somewhat with the properties 
of individual Be stars. More detailed calculations have been carried 
out (e.g. Pols et al. 1991; Raguzova 2001; de Mink et al. 2013; 
Schootemeijer et al. 2018; Hastings et al. 2021). Once more bloated 
companions are discovered, detailed comparisons with these models 
will become possible. 

Our procedure for selecting Be stars and identifying candidates for 
having a binary companion introduces some additional systematic 
uncertainty. Because the parent sample within which to search for 
Be stars began with a CMD-based selection that only includes « 
85% of Be stars (Section 2.1), some Be stars observed by APOGEE 
never enter our sample. This will not bias our results if the stars that 
are excluded have similar binary properties to the rest of the sample, 
but it could lead to an underestimate of fstripped if e.g. later-type Be 
stars or young Be stars in highly reddened regions were more likely to 
have binary companions. Be stars with sufficiently cool and luminous 
companions could also be preferentially excluded from the sample 
if these companions moved the unresolved source redward of our 
CMD search region. However, our MESA calculations predict that 
the stripped companions have low enough luminosity during most 
of the mass transfer process that this effect is minor (Figure 15). At 
present, the Poisson uncertainty on Ngetecteg = 1 is large enough to 
dominate over other sources of uncertainty. 

Our results are inconsistent at the «2 sigma level with a scenario 
in which all Be stars form via stable mass transfer in Algol-type 
systems. Not all binary evolution scenarios that could form Be stars 
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Figure 20. Constraint on the fraction of Be stars that have gone through a 
companion-stripping stage similar to HD 15124. This is inferred from the 
fact that we detected one stripped companion out of 297 Be stars and the 
timescale over which the companion is detectable with APOGEE spectra is 
z 50 times shorter than the subsequent Be star lifetime. 


go through an Algol-type phase. For example, some Be stars may 
form through binary mergers (e.g. de Mink et al. 2013; Gies et al. 
2021). This channel could produce single Be stars. 


7 CONCLUSIONS 


Motivated by suggestions that many Be stars are spun-up by mass 
transfer, we searched the spectra of Be stars from the APOGEE survey 
for stripped, slowly rotating companions. From a parent sample of 
297 Be star candidates, we identified one promising binary, HD 
15124. 

The system contains a main-sequence B star with a circumstellar 
disk (the “Be” component), and a warm, bloated subgiant companion 
(the *donor"). The subgiant nearly or completely fills its Roche lobe, 
suggesting that there is ongoing mass transfer. The Be star rotates at 
about 60% of its critical rotation velocity, faster than most normal B 
stars, but slower than most classical Be stars. We thus propose that 
its circumstellar disk is an accretion disk. Ongoing mass transfer is 
expected to continue to spin up the star and widen the orbit, such 
that the accretor will become a classical Be star when mass transfer 
ends. HD 15124 is one of a small number of well-characterized mass 
transfer binaries in this state (e.g. Harmanec et al. 2015, their Table 
7); to our knowledge, it is the first to be identified through a systematic 
search of B stars from a survey with a modelable selection function, 
enabling an estimate of occurrence rates. 

Although it satisfies the functional definition of a Be star — it is 
a B star with emission lines — we stress that the accreting star in 
HD 15124 is in a different evolutionary state from classical Be stars 
(e.g. Rivinius et al. 2013), whose emission owes to episodic mass 
loss through a viscous decretion disk. The emission lines currently 
observed in HD 15124 originate in an accretion disk, which will 
disappear once mass transfer ends. The accreting star is, however, 
being spun up by mass transfer, and our modeling assumes that when 
mass transfer ends, it will become a classical Be star. 

Modeling the system's light curves (Figures 6 and 7), high- 
resolution spectra (Figure 3, 5, 10, and 13), radial velocities (Fig- 
ure 8), spectral energy distribution (Figure 4, and 9), and astrometry 
allows us to place tight constraints on the binary's physical param- 
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eters and chemical abundances (Table 1). We interpret these using 
MESA binary evolution calculations (Figure 15, 16, 17, and 18), 
which suggest the system was born as a — (4 3) Mo binary in 
which the subgiant was initially the more massive star. Our main 
results are as follows. 


(i) Observational data: HD 15124 is a double-lined binary. The 
narrow-lined donor contributes ~ 20% of the optical flux, while the 
broad-lined Be star contributes ~ 80% (Figure 4). The circumstel- 
lar disk contributes only a few percent of the optical flux, but its 
contributions become increasingly important at longer wavelengths 
(Figure 9), and it contributes significant emission in hydrogen and 
helium lines (Figure 5). The radial velocities of the donor are well-fit 
by a circular orbit with Porb = 5.47 days and projected RV semi- 
amplitude of ~ 74km s (Figure 8). The light curve shows pho- 
tometric variability dominated by a reflection effect and ellipsoidal 
variability of the donor (Figures 6 and 7). There is also variability 
on longer and shorter timescales, likely driven by turbulence and 
changes in the structure of the disk. 

(ii) Physical parameters: We constrain the mass and radius of the 
Be star to be Mge = 5.3 + 0.6 and Rg, = 4.1 + 0.2Ro, and for the 
donor, Mgonor = 0.92 + 0.22 and Rgonor = 5.82 + 0.45 Ro (Table 1). 
The Be star is close to the main sequence (Figure 11), while the 
donor is clearly evolved. The inclination is i ~ 23 deg (close to 
face-on). The donor appears to be tidally synchronized, while the Be 
star rotates at ~60% of its critical velocity. Coupled with the fact 
that there is ongoing mass transfer, this suggests the emission lines 
originate from an accretion disk, not a decretion disk. 

(iii) Abundances: The surface abundances of both stars in HD 
15124 are quite unusual. Helium and nitrogen are enhanced by factors 
of 4 and 6, while carbon and oxygen are depleted by factors of 200 
and 6 (Figure 13). Such abundances are expected within the CNO- 
burning core of a few-solar mass star, where hydrogen has been 
converted to helium and CNO burning equilibrium is reached with 
most of the carbon converted to nitrogen. The fact that we observe 
these abundances on the stellar surface suggests that what is today 
the surface was once inside the convective core of the donor, and that 
this CNO-processed material has also been deposited on the surface 
of the Be star. 

(iv) Evolutionary history: We compared the observed properties 
of HD 15124 to binary evolution models from the BPASS library 
as well as bespoke models calculated with MESA. The models that 
most closely match HD 15124 begin with primary and secondary 
masses near 4 Mo and 3 Mo and periods of a few days. The initially 
more massive star overflows its Roche lobe shortly before or after 
the end of its main-sequence evolution, initiating a phase of thermal- 
timescale mass transfer. It evolves up the giant branch as it continues 
to lose its envelope, passing through the currently observed phase of 
HD 15124 when the mass transfer rate is around M ~ 1077 Mo yr! 
(Figures 15 and 16). In order to match the observed short orbital 
period, mass transfer is required to have been fairly non-conservative, 
with significant loss of angular momentum after Roche lobe overflow. 

Such models are predicted to terminate mass transfer within 1-2 
Myr, when almost all of the envelope has been stripped and only 
a helium core with ~ 0.1 Mo of hydrogen envelope remains. The 
donor is then predicted to contract to become a hot, compact core 
helium burning extreme horizontal branch (or “sdO/B”’) star, where 
it will be a nearly invisible companion to a rapidly rotating Be star 
spun up by accretion. This suggests that HD 15124 is a progenitor 
system for the population of Be + sdO/B binaries discovered in the 
UV (e.g. Wang et al. 2021). 

(v) Implications for the origin of Be stars: Discovery of one Be + 
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stripped star progenitor in a short-lived phase implies the existence 
of a larger population in the Be + sOB phase, when the stripped 
companion can easily escape detection. We infer that a fraction 
Stripped = Quei of Be stars have stripped companions (Fig- 
ure 20). This is broadly consistent with previous estimates from 
similar binaries containing slightly more evolved stripped stars (e.g. 
El-Badry & Quataert 2021). The fact that HD 15124 was discovered 
through a systematic search of a well-defined parent sample of Be 
stars makes this constraint more straightforwardly interpretable than 
the one inferred by El-Badry & Quataert (2021). The uncertainties 
remain large. Joint modeling of bloated and sdO/B companions to 
Be stars — many more of which will be within the reach of proposed 
future missions (e.g. Jones et al. 2021; Kulkarni et al. 2021) — will 
continue to solidify our understanding of the importance of mass 
transfer in the formation of Be stars. 
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APPENDIX A: RADIAL VELOCITIES 


Radial velocities from the APOGEE and NRES spectra are reported 
in Table A1. 


APPENDIX B: SEARCH SENSITIVITY 


Figure B1 shows noiseless simulated spectra of binaries containing 
a rapidly rotating star (representing the Be star) and a slowly rotating 
companion (the stripped star) at optical and APOGEE wavelengths. 
Our search depends on the slowly-rotating stripped star contributing 
detectable narrow lines to the spectrum at APOGEE wavelengths. 
Such companions contribute plenty of lines when the stripped com- 
panion is relatively cool, as is the case in HD 15124. However, for 
Teff, stripped 2 9,000 K, there are very few detectable metal lines in 
the H—band. The predicted composite H —band spectra are still rec- 
ognizably different from the predicted spectrum of a single rapidly- 
rotating star due to the stripped star's lower surface gravity and pro- 
jected rotation velocity. However, this difference would be challeng- 
ing to detect in real Be stars, whose hydrogen lines are contaminated 
with time-variable emission. Our search is thus sensitive primarily to 
companions with Teff, stripped S 8, 000 K; hotter companions would 
be more easily detectable in the optical. 


APPENDIX C: IS HD 15124 A TRIPLE? 


HD 15124 appears in the Hipparcos-Gaia proper motion anomaly 
catalog (Kervella et al. 2019, 2021), which consists of objects whose 
proper motion as measured by Gaia is inconsistent with the value 
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Table A1. Radial velocities of the donor 


HJD phase RV [kms7!] instrument 
2457681.8170 0.471 —25.58 + 0.96 APOGEE 
2457734.6239 0.127 —65.01 + 2.14 APOGEE 
2457762.5615 0.235 —84.82 + 2.16 APOGEE 
2457765.5750 0.786 60.56 + 1.50 APOGEE 
2458007.8932 0.091 —52.27] + 1.85 APOGEE 
2458010.8925 0.640 45.19 + 1.36 APOGEE 
2458039.8325 0.931 20.73 + 1.52 APOGEE 
2458056.7598 . 0.026 —22.95 + 0.71 APOGEE 
2458061.8121 0.950 12.03 + 1.40 APOGEE 
2458064.7898 0.494 —11.22 + 1.21 APOGEE 
2458084.6901 0.133 —67.23 + 1.62 APOGEE 
2459247.1871 0.685 56.96 + 4.90 NRES 
2459249.2283 0.058 | —45.00 + 12.85 NRES 
2459253.2535 0.794 60.52 + 5.34 NRES 
2459256.2137 . 0.335 —66.62 + 7.73 NRES 
2459257.1976 0.515 —1.80 + 6.08 NRES 
2459259.1880 0.879 38.16 + 7.40 NRES 
2459260.2069 0.065 —41.22 + 5.12 NRES 
2459260.2367 0.071 —44.74 + 4.38 NRES 
2459260.2710 0.077 -47.33 + 10.43 NRES 
2459267.1907 0.342 —62.45 + 6.41 NRES 
2459267.2539 0.354  -62.35 € 12.45 NRES 
2459267.2797 | 0.359 —66.26 + 8.38 NRES 
2459269.2240 0.714 55.56 + 6.80 NRES 
2459271.2015 0.076 —46.64 + 2.12 NRES 
2459271.2205 0.079  —54.62 + 14.74 NRES 
2459271.2351 0.082 —45.89 + 3.48 NRES 
2459272.2528 0.268 -83.40 + 14.82 NRES 
2459272.2650 | 0.270 —78.63 + 3.07 NRES 
2459274.1936 0.623 41.22 + 4.64 NRES 
2459274.2129 0.626 43.83 + 5.07 NRES 
2459274.2531 0.634 44.11 + 3.95 NRES 
2459274.2650 0.636 45.92 + 1.05 NRES 
2459275.2589 0.818 52.56 + 3.45 NRES 
2459276.2598 0.001 —13.21 + 4.20 NRES 
2459279.2124 0.540 9.59 + 2.41 NRES 
2459279.2617 0.549 13.05 + 4.13 NRES 
2459280.1959 0.720 58.42 + 2.89 NRES 
2459280.2221 0.725 58.93 + 12.30 NRES 
2459281.2054 0.905 28.50 + 5.17 NRES 
2459282.2253 0.091 —49.22 + 8.87 NRES 


inferred from the difference in the Hipparcos and Gaia positions. 
For HD 15124, the discrepancy between the Hipparcos-Gaia eDR3 
position change and the Gaia eDR3 proper motion is formally sig- 
nificant at the 15c level, with a tangential velocity anomaly of 
Av, = 1.99 + 0.13 kms^! . Some caution is warranted in interpret- 
ing this number, since neither the Hipparcos nor the Gaia astrometry 
modeled HD 15124 as an unresolved binary. However, the object's 
5-day orbital period is short compared to the Gaia and Hipparcos ob- 
servation windows, and the projected size of the binary's orbit in our 
solution is small (~ 0.2 mas) compared to the parallax (1.64 mas). 
We therefore see no strong reason to be suspicious of the reported 
proper motion anomaly. 

Assuming the proper motion anomaly is real and due to an unre- 
solved tertiary (following the methods described in El-Badry et al. 
(2021c), we do not find any resolved companions in Gaia eDR3), 
the implied mass of the unseen companion depends on its separa- 
tion. This dependence is nonlinear, because orbit-smearing effects 
are significant at periods shorter than the 25-year Hipparcos-Gaia 
epoch separation, while at periods much longer than this baseline, 
only a small fraction of the orbit is covered by the observations. For 
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Figure B1. Simulated spectra of Be + stripped star binaries at optical (left) and APOGEE (right) wavelengths. In all panels, we assume a main-sequence 
Be star with Tg, ge = 18, 000K, Rpe = 4 Ro, and v sinig, = 300km s-!. Emission lines are not included. The stripped companion has Rsuipped = 5 Ro, 
log 2stripped = 2-5, v sin istrippea = 20km s^!, and Teff, stripped increasing from top to bottom. At APOGEE wavelengths, the narrow lines that reveal the stripped 
companion's presence disappear to at Teff, stripped 2 8000 K. In this idealized case, the Brackett lines have narrower cores when a hot stripped companion is 


present, but this is unlikely to be easily recognizable in practice because these will be contaminated with time-variable emission lines from the Be star. Our 
APOGEE search is thus only sensitive to companions with Tẹ < 8000 K. Hotter companions would be detectable in the optical. 


orbital separations of 3, 5, 10, and 30 au, Kervella et al. (2021) esti- 
mate companion masses for HD 15124 of 1.6, 0.9, 0.6, and 1.0 Mo. 
Assuming a tertiary on the main sequence, a companion with any of 
these masses would be much fainter than the ~ 1000Lo Be star and 
~ 100L o donor and would easily escape detection in the spectra and 
light curves. 


Itis also worth considering again whether the Be star could be the 
distant tertiary to an inner binary containing the donor and an unseen 
companion, as is the case in v Gem (Klement et al. 2021b) and NGC 
2004 #115 (El-Badry et al. 2021a). The primary problem with this 
scenario is that it fails to explain the light curve, which shows clear 
modulation on the binary's 5.47 day orbital period (Figure 7). The 
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shape of this modulation strongly suggests it is due to a reflection 
effect: the variability is too smooth to be due to eclipses, and cannot 
be primarily ellipsoidal variation (or its period would be half the RV 
period). If it is indeed a reflection effect, the companion must be 
hotter than the donor. Accounting for dilution from the Be star, we 
see no way to hide an inner companion bright enough to cause such a 
reflection effect. Finally, we find helium enhancement in the Be star, 
which is naturally explained if it has just accreted CNO-processed 
material from the donor, but would otherwise be unusual. Although 
the existence of a third object in HD 15124 may seem contrived, it 
would not be that uncommon. Indeed, the prototypical Algol-type 
binary — Algol itself — contains a low-mass tertiary, and almost 80% 
of close binaries with Porb < 7 days have an outer tertiary (Tokovinin 
et al. 2006). About 32% of all stars in the Hipparcos catalog have a 
significant proper motion anomaly. 

We thus conclude that a faint tertiary is plausible, though the proper 
motion anomaly could also be due to systemics in the astrometry. The 
presence of a 0.5 — 2 Mo companion at a few-AU separation would 
not significantly change any aspect of our interpretation of the inner 
binary. 


This paper has been typeset from a TEX/IATEX file prepared by the author. 


MNRAS 000, 1—30 (2022) 


